Integrodifference master equation describing actively growing blood vessels in angiogenesis
aa r X i v : . [ q - b i o . T O ] S e p Integrodifference master equation describing actively growing blood vessels inangiogenesis
L. L. Bonilla, M. Carretero, and F. Terragni G. Mill´an Institute, Department of Mathematics,Universidad Carlos III de Madrid, Avenida de la Universidad 30,Legan´es 28911, Spain, e-mail: [email protected] G. Mill´an Institute, Department of Mathematics,Universidad Carlos III de Madrid, Avenida de la Universidad 30,Legan´es 28911, Spain, e-mail: [email protected] G. Mill´an Institute, Department of Mathematics,Universidad Carlos III de Madrid, Avenida de la Universidad 30,Legan´es 28911, Spain, e-mail: fi[email protected]
We study a system of particles in a two-dimensional geometry that move according to a reinforcedrandom walk with transition probabilities dependent on the solutions of reaction-diffusion equationsfor the underlying fields. A birth process and a history-dependent killing process are also considered.This system models tumor-induced angiogenesis, the process of formation of blood vessels inducedby a growth factor released by a tumor. Particles represent vessel tip cells, whose trajectoriesconstitute the growing vessel network. New vessels appear and may fuse with existing ones duringtheir evolution. Thus, the system is described by tracking the density of active tips, calculated asan ensemble average over many realizations of the stochastic process. Such density satisfies a noveldiscrete master equation with source and sink terms. The sink term is proportional to a space-dependent and suitably fitted killing coefficient. Results are illustrated studying two influentialangiogenesis models.
I. INTRODUCTION
In the present paper, we investigate the density of a set of particles undergoing a reinforced random walk, branchingand killing processes in two space dimensions. Branching and killing are point processes but the latter occurs onlywhen a particle intersects the trajectory of another particle or collides with it. These systems have received a lot ofattention as they model angiogenesis, the growth of blood vessel networks that are important in organ development,wound healing, and in many pathologies including cancer [2, 26]. However, most of the theoretical work consists ofcomputational models or mathematical models that are solved numerically and directly compared with experiments,with little mathematical elaboration [22, 27]. The goal of the present paper is to derive a master equation for thedensity of the active particles representing the angiogenic network and validate it by comparing its solutions to thoseobtained by direct simulation of the underlying stochastic process.The growth of blood vessels out of a primary vessel or angiogenesis is a complex multiscale process responsiblefor organ growth and regeneration, tissue repair, wound healing, and many other natural operations in living beings[11–13, 19, 21]. Angiogenesis is triggered by lack of oxygen (hypoxia) experienced by cells in some tissue. Such cellssecrete growth factors that diffuse and reach a nearby primary blood vessel. In response, the vessel wall opens andissues endothelial cells that move towards the hypoxic region. The cells at the tips of the advancing network arehighly motile, do not proliferate, and react to gradients of the growth factors released at the hypoxic region and tomechanical cues. Following them, other endothelial cells proliferate and build the blood vessels. Thus, capillariesbring blood, oxygen, and nutrients to the hypoxic region. Many models assume that the angiogenic network is madeout of the moving vessel tips (considered as point particles) and their trajectories. New particles are created by abranching process, while particles that intersect the trajectory of another particle are destroyed, as the vessels theyrepresent merge with a preexisting blood vessel, in a process called anastomosis. Once blood and oxygen have reachedthe hypoxic region, secretion of growth factors stops, anti-angiogenic substances may be released and a regular vesselnetwork is put in place, after pruning capillaries with insufficient blood flow.In healthy circumstances, angiogenic and anti-angiogenic activities balance. Imbalance may result in many diseasesincluding cancer [11, 17, 18]. In fact, after a tumor installed in a tissue reaches some two millimeters size, it needsadditional nutrients and oxygen to continue growing. Its hypoxic cells secrete growth factors that induce angiogenesis.Unlike normal cells, cancerous ones continue issuing growth factors and attracting blood vessels, which also supplythem with a handy transportation system to reach other organs in the body. In angiogenesis, events happeningat cellular and subcellular scales unchain endothelial cell motion and proliferation, and build millimeter size bloodsprouts and networks thereof [12, 13, 19, 21]. Models range from very simple to extraordinarily complex and often tryto illuminate some particular mechanism; see the reviews [22, 27]. Realistic microscopic models involve postulatingmechanisms and a large number of parameters that cannot be directly estimated from experiments, but they oftenyield qualitative predictions that can be tested. An important challenge is to extract mesoscopic and macroscopicdescriptions of angiogenesis from the diverse microscopic models.Early models were based on reaction-diffusion equations (RDEs) for growth factor densities, endothelial cell den-sities, etc. [14, 15, 23] and they could not describe the growing blood vessel network. In a seminal work, Andersonand Chaplain [2] derived a reinforced random walk for growing blood vessels from the RDE describing the densityof endothelial cells. Then, they supplemented the random walk with branching and anastomosis processes. Unfor-tunately, the transition probabilities obtained from a finite difference discretization of a PDE may become negativefor some parameter ranges and therefore the resulting random walk with additional branching and killing rules canbe inconsistent. On the other hand, Plank and Sleeman introduced related reinforced random walks with explicitlynon-negative transition probabilities [26]. By connecting stochastic processes to RDEs, it is possible to ascertain howdensities representing substances important for cell adhesion, anti-angiogenic substances, and other continuum fieldsinfluence a growing angiogenic network [16]. In early reinforced random walk models [2, 16, 26, 29], solutions of theRDEs are not affected by the reinforced random walk describing the advancing network. Thus, the angiogenic net-work may be seen as postprocessing the RDEs. Models using stochastic differential equations to describe blood vesselextension (instead of random walks) couple density or flux of vessel tips to RDEs for growth factors or fibronectin[4–8, 10, 30]. Therefore, concentrations of the latter are also random variables. This is also the case of multiscalemodels of vascular tumor growth that account for the angiogenesis part of the overall process by means of reinforcedrandom walks [24].The reinforced random walk setting [28] has been widely adopted to study biological systems [25]. However, toour knowledge, literature on reinforced random walk models for angiogenesis reports only numerical simulations ofvessel networks [2, 16, 22, 24, 26]. In this work, in order to quantify the process of angiogenic network growthinduced by a tumor, we define a density of active blood vessel tips, as similarly done in models based on stochasticdifferential equations [4, 7, 10, 30]. Then, in the framework of two widely studied models [2, 26], we derive adeterministic, discrete master equation for the density of active tips calculated from ensemble averages over replicasof the underlying stochastic process. Note that a density of active tips describes an advancing angiogenic network[4, 7, 10, 30] and is quite different from a density of endothelial cells or from a probability density of extendingvessel tips [2, 16, 22, 24, 26, 29]. Here, the deterministic description consists of a discrete-time integrodifferenceequation for the mentioned density coupled to discretized RDEs for the growth factor and relevant continuum fields.The most important, novel term characterizes vessel fusion and is nonlocal in time. Its space-dependent coefficientmust be estimated by comparison with numerical simulations of the involved stochastic process. Note that seekingdeterministic equations for the density of active blood vessel tips may help developing qualitative and quantitativeanalyses of more complex models of angiogenesis, such as those reviewed in [22].The paper is organized as follows. In section II, we describe the reinforced random walk and point processes modelingangiogenesis. In section III, we derive a discrete master equation for the density of active tips and compare its solutionwith the ensemble-averaged tip density calculated from the underlying stochastic process. Finally, section IV containsour conclusions and future research directions.
II. STOCHASTIC MODEL
In order to model tumor-induced angiogenesis, we first follow [2] and consider a system of nondimensional coupledRDEs for the density of endothelial cells (ECs), n ( t, x ), a growth factor (GF) concentration, C ( t, x ), and fibronectin(FN) concentration, f ( t, x ), namely ∂n∂t = D ∆ n − ∇ · (cid:18) χ αC n ∇ C (cid:19) − ∇ · ( ρ n ∇ f ) , (1) ∂C∂t = − η n C, (2) ∂f∂t = β n − γ n f. (3)In eq. (1), the flux of ECs consists of diffusion, − D ∇ n , chemotaxis proportional to the GF gradient, χn ∇ C/ (1 + αC ),and haptotaxis proportional to the FN gradient, ρn ∇ f . Diffusion arises from ECs Brownian motion, while chemotaxisindicates that ECs move in the direction of larger gradients of the GF released by the tumor. Haptotaxis is the ECsmotion to larger adhesion gradients of FN, an adhesive macromolecule attached to the extra-cellular environment.The GF in eq. (2) is simply degraded, whereas there is an uptake and production of FN as shown in eq. (3). Otherprocesses important in different situations can be included by adding more terms and new RDEs [16]. In theseequations, time has been scaled by factor τ = L /D c , where L = 2 mm is the characteristic length of the problem and D c = 2 . · − m /s is the GF diffusion coefficient. Values of the involved dimensionless parameters are given by D = 0 . χ = 0 . α = 0 . ρ = 0 . η = 0 . β = 0 .
05, and γ = 0 .
10. Note that they have been taken from [2]except for ρ = 0 .
20, which is smaller here as to simulate a weaker haptotaxis. We consider a simple square geometryFIG. 1: Sketch of the geometry for angiogenesis from a primary vessel to a circular tumor.(Figure 1), in which a primary vessel is located at x = 0 and initially emits 5 capillaries, while a small circular tumoris centered at ( x, y ) = (1 . , .
5) and acts as the GF source. Equation (1) is subject to zero-flux boundary conditionson [0 , × [0 , ξ · (cid:18) − D ∇ n + χ αC n ∇ C + ρ n ∇ f (cid:19) = 0 (4)on the boundaries of the unit square (here, ξ is the outward unit normal vector). Note that eq. (2) lacks spatialderivatives, hence the initial GF concentration has to model the tumor, as we cannot impose a (perhaps morenatural) boundary condition on the GF flux [4]. Initial conditions for n , C , and f are n (0 , x, y ) = X j =1 e − x / . − ( y − y j ) / . , (5) C (0 , x, y ) = e − (1 . − x ) / . − (0 . − y ) / . , f (0 , x, y ) = 0 . e − x / . , (6)where y = 0 . y = 0 .
5, and y = 0 . ×
201 grid of the unit square bystep h = 0 . t = 0 . n ql,m = n ( q ∆ t, lh, mh ) can be written as n q +1 l,m = n ql,m W + n ql +1 ,m W + n ql − ,m W + n ql,m +1 W + n ql,m − W , (7)with coefficients W , W , W , W , W given by W = 1 − D ∆ th + α ∆ t χ ( C ql,m )4 h (1 + αC ql,m ) h ( C ql +1 ,m − C ql − ,m ) + ( C ql,m +1 − C ql,m − ) i − ∆ t χ ( C ql,m ) h ( C ql +1 ,m + C ql − ,m + C ql,m +1 + C ql,m − − C ql,m ) − ρ ∆ th ( f ql +1 ,m + f ql − ,m + f ql,m +1 + f ql,m − − f ql,m ) , (8) W = D ∆ th − ∆ t h h χ ( C ql,m ) ( C ql +1 ,m − C ql − ,m ) + ρ ( f ql +1 ,m − f ql − ,m ) i , (9) W = D ∆ th + ∆ t h h χ ( C ql,m ) ( C ql +1 ,m − C ql − ,m ) + ρ ( f ql +1 ,m − f ql − ,m ) i , (10) W = D ∆ th − ∆ t h h χ ( C ql,m ) ( C ql,m +1 − C ql,m − ) + ρ ( f ql,m +1 − f ql,m − ) i , (11) W = D ∆ th + ∆ t h h χ ( C ql,m ) ( C ql,m +1 − C ql,m − ) + ρ ( f ql,m +1 − f ql,m − ) i , (12)where χ ( C ) = χ/ (1 + αC ) and C ql,m = C ( q ∆ t, lh, mh ), f ql,m = f ( q ∆ t, lh, mh ). Equation (7) resembles a masterequation except that W i are not true transition probabilities because they are not normalized. Let us then definetransition probabilities for a reinforced random walk as W i = W i ˜ W , with i = 0 , . . . , , ˜ W = X j =0 W j . (13) Definition 1 (reinforced random walk) Given an initial configuration with N = 5 equally spaced tips located at thegrid nodes (0 , . , (0 , . , (0 , . , (0 , . , and (0 , . , their motion is determined by means of the transitionprobabilities in eq. (13) . At each time step, a random number U (from a uniform distribution between 0 and 1) isgenerated for each tip, which will stay at the current node if U ∈ [0 , W ] , move to the left if U ∈ ( W , P j =0 W j ] ,move to the right if U ∈ ( P j =0 W j , P j =0 W j ] , move downwards if U ∈ ( P j =0 W j , P j =0 W j ] , or move upwards if U ∈ ( P j =0 W j , . The transition probabilities are evaluated at the tip node, movement is allowed only onto the nodesof the spatial grid and all tips reaching any of the four boundaries of the unit square are ‘deactivated’. Definition 2 (anastomosis) Whenever a tip meets an existing vessel (namely, the trajectory of another tip), it mergeswith the latter and is ‘deactivated’. Whenever a tip meets another tip, one of them is ‘deactivated’ and the other oneremains active.
Definition 3 (branching) A new tip can branch out of any of the existing active tips with a certain probability,regardless of its age. This is different from [2], in which a minimum age is required for a tip to branch. Tip branchingis simulated as follows. At a given time step and for an active tip at the node located at x = ( x, y ) , a random number U b is extracted from a uniform distribution between and . . Then, a new tip branches out of that tip at the samenode if U b ≤ φ ( x ) ∆ t , where φ ( x ) = if ≤ x ≤ . . if . < x ≤ . . if . < x ≤ . if . < x ≤ . if . < x ≤ . (14) The upper bound of U b produces a nonzero probability that no tip branches out at a given time. Note that a newlybranched tip will move according to Definition 1 and that branching onto an occupied node of the grid or any of theunit square boundaries is not allowed. In the same context, Plank and Sleeman proposed a different reinforced random walk that, in absence of branchingand anastomosis, yields in the continuum limit the PDE in eq. (1) for the probability density [26]. The associatedtransition probabilities are ˆ τ H ± l,m = 4 Dh τ ( w l ± ,m ) τ ( w l + ,m ) + τ ( w l − ,m ) + τ ( w l,m + ) + τ ( w l,m − ) , (15)ˆ τ V ± l,m = 4 Dh τ ( w l,m ± ) τ ( w l + ,m ) + τ ( w l − ,m ) + τ ( w l,m + ) + τ ( w l,m − ) , (16)where w = ( C, f ) and τ ( w ) = τ ( C ) τ ( f ), with τ ( C ) = (1 + αC ) χαD , τ ( f ) = e ρf/D . (17)Note that in eqns. (15)-(16) transitions are calculated by means of the nearest half-step neighbors instead of thenearest neighbors [25]. Tip motion, anastomosis, and branching can be treated as in Definitions 1, 2, and 3 above.Numerical simulations of the stochastic model with transition probabilities given by either eq. (13) or eqns. (15)-(16)produce angiogenic networks that are quite similar. A typical outcome is given in Figure 2, which shows the growingof a vessel network on the underlying grid from the primary vessel towards the tumor. x y (a) x y (b) (c) x y (d) x y FIG. 2: Growing vessel network on the underlying grid simulated with transition probabilities as in eqns. (15)-(16)at (a) 5 days (35 active tips), (b) 6 days (36 active tips), (c) 7 days (52 active tips), (d) 8 days (38 active tips).
III. MASTER EQUATION
In the configuration described in section II, vessel tips tend to cluster in a relatively narrow region due to thegradients of the GF and FN. Thus, many of them are killed by anastomosis during the evolution and the number ofactive tips maintains rather small. As a consequence, the stochastic model is not self-averaging and average quantitiesare not expected to resemble those of a typical realization (i.e., a replica) of the process. Indeed, different replicasof the stochastic process produce angiogenic networks (i.e., sets of the trajectories of all actively moving vessel tips)that may look quite different. Nevertheless, the tip density is preserved under the ensemble average over a sufficientlylarge number of replicas and thus a deterministic description can follow for ensemble-averaged quantities. Let us thendefine the density of active tips at time t as an ensemble average over realizations of the stochastic process. Definition 4 (ensemble-averaged tip density) Consider all N ( t, ω ) active tips at time instant t for a realization ofthe stochastic process labeled by ω , with ω = 1 , . . . , N . Let X i ( t, ω ) be the node-location of the i -th tip at time instant t . Then, the ensemble-averaged tip density P N ( t, x ) is defined as P N ( t, x ) = 1 N N X ω =1 N ( t,ω ) X i =1 δ σ ( x − X i ( t, ω )) , δ σ ( x ) = e − x / ( σ x ) − y / ( σ y ) πσ x σ y . (18)Whenever it exists, the limit of P N ( t, x ) as N → ∞ and σ x , σ y → P ( t, x ). In oursimulations, we use standard deviations equal to σ x = 0 .
03 and σ y = 0 .
05, being results fairly robust in connectionwith such values. Figure 3 shows time evolution of the ensemble-averaged tip density with N = 800, when usingtransition probabilities as in eqns. (15)-(16) in the stochastic model. Note that, once P N ( t, x ) is computed, theensemble-averaged number of active tips can also be recovered by integrating the density over the unit square andtaking the integer part (this is used, for a better comparison with the master equation solution described below, after P N ( t, x ) has reached the tumor).FIG. 3: Ensemble-averaged (over N = 800 replicas) stochastic tip density, when using transition probabilities as ineqns. (15)-(16), at (a) 5 days (19 active tips), (b) 6 days (32 active tips), (c) 7 days (41 active tips), (d) 8 days (29active tips).We now derive a master equation for the density of active vessel tips corresponding to the stochastic processdefined in section II. The master equation should be discrete in time and contain the transition probabilities, fromeither eq. (13) or eqns. (15)-(16). In addition, it should include terms related to branching and anastomosis. Let P ql,m = P ( q ∆ t, lh, mh ) be the density of active tips at the grid node indexed by ( l, m ) and time t = q ∆ t . Then, thediscrete master equation associated with transition probabilities as in eq. (13) is given by P q +1 l,m = P ql,m W + P ql +1 ,m W + P ql − ,m W + P ql,m +1 W + P ql,m − W + ∆ t φ l,m P ql,m − ∆ t Γ l,m P ql,m q X k =0 P kl,m . (19)The last-but-one is a birth term modeling branching of new tips, where φ l,m is the discrete counterpart of the branchingfunction given in eq. (14). The last one is a nonlocal killing term corresponding to anastomosis, where Γ l,m is thediscrete anastomosis function to be fitted by comparison with the ensemble-averaged results of the stochastic model,as discussed below. The anastomosis coefficient function depends on the spatial grid because the branching functiondoes. Similar birth and death terms have been incorporated to Fokker-Planck-type equations in angiogenesis modelsthat describe vessel extension by stochastic differential equations instead of random walks [4, 9, 30]. It is worthremarking that the inclusion of both source terms in eq. (19) is essential in order to properly describe an advancingvessel network and further analyze its behavior [4–6, 8, 9, 30]. Remark 1
The transition probabilities W i , with i = 0 , . . . , , depend on the solutions of eqns. (1) – (3) calculated usingthe explicit Euler scheme as discussed in section II. Thus, the evolution of the active vessel tips does not affect thedensity of ECs and the concentration of GF or FN. Integrodifference eq. (19) has to be solved with absorbing boundary conditions, P ql,m = 0, at all nodes on the fourboundaries of the unit square [20]. The initial density of active tips is set to P l,m = 2 πσ x σ y e − x l / ( σ x ) N X i =1 e − ( y m − Y i ) / ( σ y ) , (20)where N = 5, Y , . . . , Y are the ordinates of the initial tips (see Definition 1), and ( x l , y m ) = ( lh, mh ). Remark 2
During time evolution, the number of active tips is set equal to N = 5 as long as the x -location of themaximum of the tip density is smaller than . (i.e., in absence of both branching and anastomosis, see eq. (14) andeq. (22) below). Otherwise, due to the structure of P l,m in eq. (20) , N is the integer part of the sum P l,m P ql,m h over all grid nodes at each time instant t = q ∆ t . In the limit h → D/h → ∞ , the solution of eq. (19) should approach the ensemble-averaged tip densityprovided by the stochastic model, namely P N for sufficiently large N (see Figure 3). Remark 3
Numerical simulations show that W in eq. (13) may become negative during some time intervals. Ac-cording to Definition 1 of the reinforced random walk, the tip cannot move to the left because W > W + W for thesetimes. Thus motion to the left is impeded in favor of motion to the right and the resulting random motion toward thetumor is artificially accelerated. By construction, W depends on the values of D , χ , α , ρ (together with the GF and FN derivatives), as seen ineqns. (8)–(13). As W turns out to be negative during various iterations of the stochastic simulations, the intervalsover which we select probabilities are inconsistently defined, which privileges advance toward the tumor. Hence,the resulting random motion will not be a true random walk, but it may still approximate well the solution of thePDE in eq. (1). Therefore, the transition rules in Definition 1 become inconsistent for some ranges of values of thementioned quantities, which artificially privileges advance of the vessels toward the tumor. For the parameter valuesconsidered in this paper (which are the same as in [2]), the discrete master equation solution cannot be matched withthe outcome of the reinforced random motion induced by eq. (13) plus branching and anastomosis.Now, in order to find a deterministic description that agrees with ensemble averages of the stochastic processdiscussed in section II, let us take into account the (always positive) transition probabilities given by eqns. (15)-(16).The resulting discrete master equation is eq. (19) with coefficients W , . . . , W given by W = 1 − (ˆ τ H + l,m + ˆ τ H − l,m + ˆ τ V + l,m + ˆ τ V − l,m )∆ t, W = ˆ τ H − l +1 ,m ∆ t, W = ˆ τ H + l − ,m ∆ t, W = ˆ τ V − l,m +1 ∆ t, W = ˆ τ V + l,m − ∆ t, (21)and initial condition as in eq. (20). Note that Remarks 1 and 2 above still hold for the new discrete master equation (19)based on eq. (21). Its resolution can be performed using the same spatial grid and ∆ t = 0 .
001 as in the correspondingstochastic process. Similarly to the work carried out in models based on stochastic differential equations [30], thediscrete anastomosis function Γ l,m should be estimated so that the number of active tips calculated by ensembleaverages of the stochastic model and by numerically solving the discrete master equation agree. Here, Γ l,m is thediscrete counterpart of the function Γ( x ) = ≤ x ≤ . . . < x ≤ . . . < x ≤ . . . < x ≤ . . . < x ≤ , (22)which has been fitted in such a way that the maximum difference between the total number of active tips calculatedvia eqns. (19)–(21) and by stochastic ensemble averages based on eqns. (15)-(16) is equal to two for each time instant. t (days) t i p s nu m b e r FIG. 4: Time evolution of the total number of active tips as computed by the ensemble-averaged (over N = 800replicas) stochastic process based on eqns. (15)-(16) (blue line) and the discrete master equation (19) with eq. (21)(red line).Figure 4 shows the time evolution of the total number of active tips according to the solution of the master equationand ensemble averages of stochastic simulations based on eqns. (15)-(16). Both descriptions agree rather well. Thenumber of active tips reaches a maximum and then decrease after the first ones arrive at the tumor.FIG. 5: Tip density from the discrete master equation (19) with eq. (21) at (a) 5 days (20 active tips), (b) 6 days(33 active tips), (c) 7 days (41 active tips), (d) 8 days (31 active tips).We have solved the discrete master equation (19) with transition probabilities given by eq. (21) for the anastomosisfunction given by eq. (22). The resulting density of active tips evolves in time as depicted in Figure 5. This figureshould be compared to the evolution of the active tip density calculated from ensemble averages and exhibited inFigure 3. Note that the advance of the moving lump of active vessel tips is quite similar in both figures, althoughthe shapes and sizes of the moving lumps are not the same. At the end of the formation stage, the lump obtained byensemble averages of the random walk model is quite similar to that obtained by solving the master equation (19),as shown in Figures 3(a) and 5(a). Later on, the lump that solves (19), shown in Figures 5(b)-(d), becomes narrowerthan that in Figures 3(b)-(d) while still keeping the same overall number of tips, as depicted in Figure 4. Therefore,the corresponding density inside the lump has to be larger for the solution of the master equation in Figure 5 than forthe ensemble averages of trajectories shown in Figure 3: the maximum density for the solution of the master equationmay become almost twice the ensemble-averaged maximum density. t (days) t i p s nu m b e r FIG. 6: Comparison of the total number of active tips calculated from ensemble averages (over N = 800 replicas) ofthe stochastic process based on eq. (13) (red curve) and its counterpart based on eqns. (15)-(16) (blue curve).Figure 6 compares the ensemble-averaged total number of active tips when using the transition probabilities givenby either eq. (13) or eqns. (15)-(16). This number decreases as the advancing network reaches the tumor. Thus,Figure 6 highlights that tips moving according to eq. (13) arrive at the tumor before (about 0 . W in eq. (13) becomes negative during the evolution of the angiogenic network. As explained in Remark3, W < ≈ √ τ = τ / √
2; see section II).The continuum limit of the discrete integrodifference master equation (19) based on eq. (21) is ∂P∂t + ∇· (cid:18) χ ∇ C αC P (cid:19) + ∇· ( ρ ∇ f P ) − D ∆ P = φ ( x ) P − Γ( x ) P Z t P ( s, x ) ds. (23)This integrodifferential equation is similar to that for the marginal density of active tips derived in [5, 6] for a lattice-free model based on Langevin-Ito stochastic differential equations [4, 30] instead of reinforced random walks. Anumerical finite difference scheme to solve eq. (23) is studied in detail in [9]. There it is proved that the solutions ofthe scheme are positive, stable and converge to the solutions of the integrodifferential equation. The main differencesbetween eq. (23) and the equations analyzed in those works are that φ ( x ) was a function of the GF concentration andΓ( x ) was a constant. These differences stem from the different definition of branching used in the model of [5, 6, 9, 30](based on Langevin-Ito stochastic differential equations) and in the model considered here (based on [2, 26]). Onthe other hand, numerical simulations of the stochastic process and solutions of the discrete master equation showthat the density of active tips forms a moving lump after an initial stage. The profile of such lump at y = 0 . IV. CONCLUSION AND OUTLOOK
In this paper, we have derived for the first time a discrete master equation for the density of active blood vesseltips starting from two well-known models of angiogenesis [2, 26]. These and other models based on reinforced randomwalks [16, 24] postulate that the densities of endothelial cells, growth factor, and appropriate substances solve asystem of coupled reaction-diffusion equations. They propose reinforced random walks (either occurring on a latticeor being lattice-independent) whose continuum limit is the partial differential equation for the density of endothelialcells. Then, tip branching and vessel fusion are added to the random walks. We have shown that the evolutionof the angiogenic network is described by the ensemble-averaged density of active tips. The latter solves a discretemaster equation with birth (branching) and killing (anastomosis) terms, which are absent in all previous works. It0is important to ensure that the terms corresponding to transition probabilities in this equation are non-negative, asdone in [26]. Otherwise, the solution of a master equation obtained from any numerical finite difference scheme (asdone in [2]) may evolve differently from ensemble averages of the active tips. In the continuum limit, this masterequation becomes an integrodifferential equation similar to that describing the density of active tips in Langevin-Itomodels [4, 5, 30]. Rigorously establishing our derivations, and proving existence, uniqueness, and long-time stabilityresults for the limiting master equation would be highly desirable, but clearly out of the scope of the present paper.There are different interesting directions in which this work can be extended. The density of active tips satisfiesan equation that is alike that of the density of endothelial cells except for the addition of source and sink terms. Thelatter is nonlocal in time. One interesting direction is obtaining a soliton description for the moving lump of Figures 3and 5, as similarly done in [5, 6, 8]. Finding equations for the fluctuations about the density of active tips (whichis an average quantity) would be worthwhile also. The fluctuations would theoretically describe the statistics of theangiogenic network, a topic which is quite unexplored at the present time.Another direction is extending the methodology of this paper to models of angiogenesis that probe cellular lengthscales. For example, cellular Potts models directly describe cells that change size and the extracellular matrix byMonte Carlo dynamics coupled to continuum fields for growth factors concentration, elastic fields, etc. (see reviews[22, 27] and articles [3, 31, 32]). It is challenging to deduce dynamics for the active vessel tips from the Monte Carlocellular dynamics. Once this is done, it should be possible to derive a master equation with source terms for the activetip density. This could be more productive than deriving an equation for the density of endothelial cells with a givencenter of mass and a given perimeter, as in [1].
ACKNOWLEDGMENTS
This work has been supported by the FEDER/Ministerio de Ciencia, Innovaci´on y Universidades – Agencia Estatalde Investigaci´on grant MTM2017-84446-C2-2-R. [1] M. Alber, N. Chen, P. M. Lushnikov and S. A. Newman, Continuous macroscopic limit of a discrete stochastic model forinteraction of living cells, Phys. Rev. Lett. (2007), 168102.[2] A. R. A. Anderson and M. A. J. Chaplain, Continuous and discrete mathematical models of tumor-induced angiogenesis,B. Math. Biol. (1998), 857–900.[3] A. L. Bauer, T. L. Jackson and Y. Jiang, A cell-based model exhibiting branching and anastomosis during tumor-inducedangiogenesis, Biophys. J. (2007), 3105–3121.[4] L. L. Bonilla, V. Capasso, M. Alvaro and M. Carretero, Hybrid modeling of tumor-induced angiogenesis, Phys. Rev. E (2014), 062716.[5] L. L. Bonilla, M. Carretero, F. Terragni and B. Birnir, Soliton driven angiogenesis, Sci. Rep. (2016), 31296.[6] L. L. Bonilla, M. Carretero and F. Terragni, Solitonlike attractor for blood vessel tip density in angiogenesis, Phys. Rev.E (2016), 062415.[7] L. L. Bonilla, V. Capasso, M. Alvaro, M. Carretero and F. Terragni, On the mathematical modelling of tumor-inducedangiogenesis, Math. Biosci. Eng. (2017), 45–66.[8] L. L. Bonilla, M. Carretero and F. Terragni, Ensemble averages, soliton dynamics and influence of haptotaxis in a modelof tumor-induced angiogenesis, Entropy (2017), 209.[9] L. L. Bonilla, A. Carpio, M. Carretero, G. Duro, M. Negreanu and F. Terragni, A convergent numerical scheme forintegrodifferential kinetic models of angiogenesis, J. Comput. Phys. (2018), 1270–1294.[10] V. Capasso and D. Morale, Stochastic modelling of tumour-induced angiogenesis, J. Math. Biol. (2009), 219–233.[11] P. Carmeliet, Angiogenesis in life, disease and medicine, Nature (2005), 932–936.[12] P. Carmeliet, M. Tessier-Lavigne, Common mechanisms of nerve and blood vessel wiring, Nature (2005), 193–200.[13] P. Carmeliet and R. K. Jain, Molecular mechanisms and clinical applications of angiogenesis, Nature (2011), 298–307.[14] M. A. J. Chaplain and A. Stuart, A model mechanism for the chemotactic response of endothelial cells to tumour angio-genesis factor, IMA J. Math. Appl. Med. (1993), 149–168.[15] M. A. J. Chaplain, The mathematical modelling of tumour angiogenesis and invasion, Acta Biotheor. (1995), 387–402.[16] M. A. J. Chaplain, S. R. McDougall and A. R. A. Anderson, Mathematical modeling of tumor-induced angiogenesis, Annu.Rev. Biomed. Eng. (2006), 233–257.[17] J. Folkman, Tumor angiogenesis: therapeutic implications, New Engl. J. Med. (1971), 1182–1186.[18] J. Folkman, Angiogenesis, Annu. Rev. Med. (2006), 1–18.[19] M. Fruttiger, Development of the retinal vasculature, Angiogenesis (2007), 77–88.[20] C. W. Gardiner, Stochastic methods. A handbook for the natural and social sciences, Springer-Verlag Berlin, Heidelberg,2010.[21] R. F. Gariano and T. W. Gardner, Retinal angiogenesis in development and disease, Nature (2005), 960–966. [22] T. Heck, M. M. Vaeyens and H. Van Oosterwyck, Computational models of sprouting angiogenesis and cell migration:towards multiscale mechanochemical models of angiogenesis, Math. Model. Nat. Pheno. (2015), 108–141.[23] L. A. Liotta, G. M. Saidel and J. Kleinerman, Diffusion model of tumor vascularization, B. Math. Biol. (1977), 117–128.[24] P. Macklin, S. McDougall, A. R. A. Anderson, M. A. J. Chaplain, V. Cristini and J. S. Lowengrub, Multiscale modellingand nonlinear simulation of vascular tumour growth, J. Math. Biol. (2009), 765–798.[25] H. G. Othmer and A. Stevens, Aggregation, blowup and collapse: the ABC’s of taxis and reinforced random walks, SIAMJ. Appl. Math. (1997), 1044–1081.[26] M. J. Plank and B. D. Sleeman, Lattice and non-lattice models of tumour angiogenesis, B. Math. Biol. (2004), 1785–1819.[27] M. Scianna, J. Bell and L. Preziosi, A review of mathematical models for the formation of vascular networks, J. Theor.Biol. (2013), 174–209.[28] F. Spitzer, Principles of random walk, Springer-Verlag New York, New York, 2001.[29] A. St´ephanou, S. R. McDougall, A. R. A Anderson and M. A. J. Chaplain, Mathematical modelling of the influence ofblood rheological properties upon adaptative tumour-induced angiogenesis, Math. Comput. Model. (2006), 96–123.[30] F. Terragni, M. Carretero, V. Capasso and L. L. Bonilla, Stochastic model of tumor-induced angiogenesis: ensembleaverages and deterministic equations, Phys. Rev. E (2016), 022413.[31] R. F. M. Van Oers, E. G. Rens, D. J. La Valley, C. A. Reinhart-King and R. M. H. Merks, Mechanical cell-matrix feedbackexplains pairwise and collective endothelial cell behavior in vitro, PLoS Comput. Biol. (2014), el003774.[32] R. Vega, M. Carretero, R. D. M. Travasso, and L. L. Bonilla, Notch signaling and taxis mechanims regulate early stageangiogenesis: A mathematical and computational model, PLoS Comput. Biol.16