Comparative analysis of continuum angiogenesis models

Although discrete approaches are increasingly employed to model biological phenomena, it remains unclear how complex, population-level behaviours in such frameworks arise from the rules used to represent interactions between individuals. Discrete-to-continuum approaches, which are used to derive systems of coarse-grained equations describing the mean-field dynamics of a microscopic model, can provide insight into such emergent behaviour. Coarse-grained models often contain nonlinear terms that depend on the microscopic rules of the discrete framework, however, and such nonlinearities can make a model difficult to mathematically analyse. By contrast, models developed using phenomenological approaches are typically easier to investigate but have a more obscure connection to the underlying microscopic system. To our knowledge, there has been little work done to compare solutions of phenomenological and coarse-grained models. Here we address this problem in the context of angiogenesis (the creation of new blood vessels from existing vasculature). We compare asymptotic solutions of a classical, phenomenological “snail-trail” model for angiogenesis to solutions of a nonlinear system of partial differential equations (PDEs) derived via a systematic coarse-graining procedure (Pillay et al. in Phys Rev E 95(1):012410, 2017. https://doi.org/10.1103/PhysRevE.95.012410). For distinguished parameter regimes corresponding to chemotaxis-dominated cell movement and low branching rates, both continuum models reduce at leading order to identical PDEs within the domain interior. Numerical and analytical results confirm that pointwise differences between solutions to the two continuum models are small if these conditions hold, and demonstrate how perturbation methods can be used to determine when a phenomenological model provides a good approximation to a more detailed coarse-grained system for the same biological process. Supplementary Information The online version contains supplementary material available at 10.1007/s00285-021-01570-w.


Introduction
Continuum and discrete approaches are widely used to develop theoretical models for biological systems. In the former framework, one approximates quantities of interest using continuous variables and measures their evolution in space and/or time via differential equations. By contrast, discrete approaches describe the behaviour of individual "agents" (such as cells or molecules) in simulations Van Liedekerke et al. 2015;Metzcar et al. 2019). These two modelling frameworks are not mutually exclusive: for instance, "hybrid models" use continuous approaches to represent certain quantities of interest (e.g., the concentration of a drug), but apply discrete methods for other quantities such as cells (Rejniak and Anderson 2011). Hybrid and discrete models for a biological process can typically incorporate more faithful representations of underlying cellular and/or molecular mechanisms than phenomenological continuum systems, and they are becoming increasingly prevalent in the mathematical biology community (Van Liedekerke et al. 2015;Metzcar et al. 2019).
Despite the advantages and wide use of agent-based models (ABMs), we still have a poor understanding of how different parameter values and agent interaction rules affect their solutions over large spatial and temporal scales. Kursawe et al. (2017), for example, found that results from a particular discrete model (the vertex model) are sensitive to the manner in which it is numerically implemented. Another study by Osborne et al. (2017) compared discrete approaches for modelling cell proliferation, adhesion, and signalling using a single computational framework. The models that they investigated (which included cellular automata, the Cellular Potts model, the overlapping spheres model, Voronoi tessellations, and the vertex model) did not always generate equivalent summary statistics; their work thereby illustrated how inferred biological conclusions may depend on the modelling approach used to generate results. Although such studies represent a much needed first step for evaluating discrete approaches and their underlying rules, it is challenging to gain rigorous insight into the general behaviour of ABMs unless there is a common mathematical framework in which to analyse them.
Coarse-grained macroscopic partial differential equations (PDEs) represent a possible candidate for this framework. They are constructed from discrete-to-continuum derivations, rather than from phenomenological arguments, and describe how the average distribution of agents in a discrete model evolves over time and space (Othmer et al. 1988;Baker et al. 2010;Simpson and Baker 2011;Markham et al. 2013;Bonilla et al. 2016;Dyson and Baker 2015;Spill et al. 2015;Buttenschön et al. 2018;Motsch and Peurichard 2018;Chaplain et al. 2020). Since continuum models are often easier to analyse and simulate than discrete ones, they can be used to (indirectly) determine how the microscopic rules of a discrete model generate complex collective dynamics. However, mean-field models derived from discrete-to-continuum approaches tend to be highly nonlinear, and the degree of nonlinearity depends on the rules in the underlying discrete model (Simpson et al. 2009;Penington et al. 2011;Bruna and Chapman 2012;Penington et al. 2014;Pillay et al. 2018). To our knowledge, little recent attention has focussed on the systematic comparison of nonlinear mean-field models derived from discrete-to-continuum approaches (Horstmann et al. 2004). In this article, we aim to shed light on the connection between coarse-grained and phenomenological models by analysing solutions from two particular continuum systems for angiogenesis (the growth of new blood vessels from pre-existing vasculature). Angiogenesis is implicated in many essential biological processes, such as wound healing and tissue growth, because the resulting blood vessels provide a source of nutrients for damaged and/or developing tissue (Potente et al. 2011;Viallard and Larrivée 2017;Duran et al. 2017; Barrientos et al. 2008;Risau 1997). Angiogenesis also constitutes a crucial step in cancer development, since it supplies cancerous cells with nutrients for continued growth, creates a system to remove metabolic waste, and facilitates metastasis, which is the invasion of tumour cells into new regions of tissue (Folkman 1995;Byrne 2010;Carmeliet and Jain 2011;Duran et al. 2017;Viallard and Larrivée 2017).
The models we consider in this article describe the dynamics of two cell types involved in angiogenesis: tip and stalk cells (Duran et al. 2017;Betz et al. 2016;Viallard and Larrivée 2017;Carmeliet and Jain 2011;Phng and Gerhardt 2009;Gerhardt et al. 2003;Hellström et al. 2007). Tip cells guide the migration of new sprouts by moving chemotactically up gradients of angiogenic factors (also known as tumour angiogenic factors, or TAFs), which include such molecules as vascular endothelial growth factor (Yadav et al. 2015;Lin et al. 2016). Aside from migrating through the extracellular matrix, tip cells can also branch to create multiple angiogenic sprouts and fuse with other tip or stalk cells, in a process called anastomosis, to create closed loops in the network. Stalk cells, meanwhile, have a more proliferative phenotype than tip cells and are located along the path of tip cell migration. They provide structural support to the new network, and help to establish a lumen through which blood can flow (Duran et al. 2017;Betz et al. 2016;Szymborska and Gerhardt 2018;Iruela-Arispe and Davis 2009).
A variety of continuum and discrete approaches have been used to model tip and stalk cell dynamics; they include (but are not limited to) compartment-based models (Spill et al. 2015), cellular automata (Anderson and Chaplain 1998;Qutub and Popel 2009;Jackson and Zheng 2010;Pillay et al. 2017), the Cellular Potts model (Bauer et al. 2007;Boas et al. 2018), deterministic PDEs (Flegg et al. 2015;Connor et al. 2015), and stochastic differential equations (Bonilla et al. 2014(Bonilla et al. , 2016Perfahl et al. 2017;Bonilla et al. 2020). They may be implemented in specialised computational frameworks such as Microvessel Chaste or CompuCell3D, to name but two (Grogan et al. 2017;Swat et al. 2009). For a more detailed list of mathematical models for angiogenesis, we refer the interested reader to the following reviews (Scianna et al. 2013;Mantzaris et al. 2004;Heck et al. 2014;Flegg et al. 2020).
The first continuum model that we consider is a coarse-grained system derived from a rule-based ABM in Pillay et al. (2017). The Pillay et al. ABM (henceforth denoted P-ABM) tracks the movement and proliferation of tip and stalk cells in response to a generic TAF on a 2D unit square lattice. Tip cells move via a biased random walk towards increasing TAF concentrations and proliferate to create new sprout branches.
Stalk cells are created at sites left vacant when tip cells move, and anastomosis occurs when a tip cell moves into an occupied lattice site.
The coarse-grained PDE describes the P-ABM's average distribution of tip and stalk cells over time and space. Although the full system of equations is formulated in 2D, it may be reduced to one spatial variable by column averaging the dependent variables, making a mean-field approximation, and non-dimensionalising the equations with appropriate length and time scales. The resulting 1D coarse-grained system can be written in Cartesian coordinates as follows: we denote by N (x, t) and E(x, t) the column averaged tip and stalk cell densities, respectively, at location x ∈ [0, L] and time t ∈ [t 0 , ∞), with t 0 ≥ 0, and by C(x, t) the column averaged concentration of a generic TAF that regulates the cell dynamics. The system is closed by imposing no-flux boundary conditions for the tip cells (boundary conditions are not required for the stalk cell equation because it is a first-order differential equation with respect to t -see below). We assume that the source of TAF is located at x = L > 0, that there is no TAF at the original vessel where x = 0, and that the decay of TAF and its uptake by cells are negligible. Under these assumptions, the coarse-grained model is given by production from tip cell movement + a n μN 2 + a n N D where D > 0 is the tip cell random motility coefficient, χ > 0 represents tip cell sensitivity to TAF, λ > 0 is related to the rate of tip cell branching, μ > 0 is related to tip cell movement, and the non-negative functions G(x) and H (x) represent the tip and stalk cell densities at time t 0 ≥ 0, respectively. The values of D, χ , μ, and λ may be written in terms of the P-ABM parameters (Pillay et al. 2017). In Eq. (1), the parameter 0 ≤ a e ≤ 1 is related to the prevalence of tip-to-sprout anastomosis, and in practice is estimated by fitting to data from the P-ABM using a nonlinear least squares method. The parameter a n , by contrast, is a model choice parameter: it either takes the value 0 or 1, depending on whether tip-to-tip anastomosis is included in the microscopic model. For the remainder of this article, we fix a n = 1. We refer to Eqs.
(1)-(4) as the Pillay PDE system, or P-PDE for short (see Table 1). In Eq. (1), the expressions representing tip cell movement are multiplied by the nonlinear term (1 − a n N − a e E), which derives from the volume exclusion rules incorporated in the P-ABM and models the effects of cell elimination due to tipto-tip and tip-to-sprout anastomosis. This term, however, may cause the P-PDE to become ill-posed if a n N + a e E > 1 (Pillay et al. 2017(Pillay et al. , 2018. When this occurs, the coarse-grained tip cell equation becomes similar to a negative diffusion model, which is known to yield ill-posed solutions that do not depend continuously on the initial data. Thus, the PDE solutions have the potential to be ill-posed if there are sufficiently large densities of tip and stalk cells. Although in practice multiple signalling cues act to prevent the emergence of such high cell density regions, for instance by creating a "salt-and-pepper" pattern of tip cells as they emerge from a pre-existing vessel (Bentley et al. 2008;Blanco and Gerhardt 2013;Duran et al. 2017), these subcellular processes are outside the scope of this coarse-grained model. Therefore, to ensure well-posedness of the PDE solutions, we simulate the system from t 0 > 0, by which time regions of high tip cell density have been sufficiently eliminated by anastomosis so that (1 − a n N − a e E) ≥ 0. Guided by previous numerical implementations from Pillay et al. (2017), we choose t 0 = 0.2 with column averaged P-ABM results as initial conditions (see Appendix E for details on the numerical methods).
Equation (3) describes the TAF dynamics. Byrne and Chaplain (1995) determined that for typical parameter values 0 < γ 1, hence a quasi-steady state approximation can be made for C(x, t) by setting γ = 0. The solution to Eq. (3) in this limit, subject to boundary conditions (4), is C(x, t) ≈ C(x) = νx, where ν = L −1 . We will assume for the remainder of the article that the TAF concentration is given by this expression.
The second continuum angiogenesis model that we investigate is based on the classical "snail-trail" framework. Its defining feature is that the stalk cell density increases at a rate proportional to the magnitude of net tip cell flux, which models the experimental observation that stalk cells proliferate along the path of tip cells (Balding and McElwain 1985;Byrne and Chaplain 1995;Connor et al. 2015). Although originally constructed from a phenomenological scheme proposed by Balding and McElwain (1985), snail-trail models may also be derived from discrete, compartment-based models for angiogenesis (Spill et al. 2015). In these PDE systems, tip cells are assumed to move via chemotaxis up TAF gradients and by random motion. Other processes can also be modelled, such as branching of new tip cells and tip cell elimination by anastomosis (Byrne and Chaplain 1995;Pettet et al. 1996). If N (x, t), E(x, t), C(x, t), D, χ , a n , a e , μ, and λ have the same physical interpretations as in the P-PDE, then the non-dimensional 1D version of the snail-trail model may be written in Cartesian coordinates as production from net tip cell flux (6)

Table 1
The two continuum angiogenesis models investigated in this article

Name of model
Tip and stalk cell evolution equations (stalk cell density), C(x, t) (TAF concentration, governed by Eq. (3), whose solution is C(x, t) The biological significance of each parameter value is described in the text. All models are closed with the initial and boundary conditions given by Eq. (4). The function The initial conditions, boundary conditions, and TAF dynamics for this system are given by Eqs.
(3)-(4). We again make a quasi-steady state approximation for the TAF dynamics to conclude that the TAF concentration in the snail-trail model is given by The absolute value sign in Eq. (6) ensures that the stalk cell density is always nonnegative, even when the net tip cell flux is negative. In Martinson et al. (2020), we argued that the factor κ(x) must be included in the snail-trail framework as a correction for using a net quantity (namely, the net tip cell flux) to estimate the total stalk cell density rate of change. Its derivation, which is reproduced in Appendix F, makes use of the P-ABM rules to obtain the following approximate formula for κ(x), which is valid in situations where the TAF gradient is non-zero: The corrective factor is a function of x because the TAF gradient may vary with respect to space. We have already noted, however, that in this article our column averaged TAF field is C(x) = νx; since the gradient of this field is constant, this means that κ(x) = μ/(χ ν) is also constant. We will refer to Eqs. (5)-(6) as the snail-trail model (ST-PDE, see Table 1). Although the P-PDE has a different stalk cell evolution equation than the ST-PDE and contains additional nonlinear terms, Pillay et al. (2017) found that, for certain parameter regimes, numerical solutions to the two continuum models were indistinguishable. We observe this behaviour in Fig. 1, where we present numerical solutions to the ST-PDE and P-PDE. In Fig. 1a, we find that at each time point shown, the largest pointwise difference between solutions is less than 5% of the maximum tip cell density. Similarly, in Fig. 1b the largest pointwise difference between stalk cell solutions is less than 3% of the maximum stalk cell density for each time point shown. Results such as these raise the question of whether solutions to the ST-PDE and P-PDE agree for larger regions of parameter space. Addressing this issue is important, as it would establish conditions under which nonlinear coarse-grained models such as the P-PDE can be represented well by simpler phenomenological systems like the ST-PDE.

Article outline
In this article, we compare the behaviour of solutions to the ST-PDE and P-PDE in order to determine when they will be in good agreement. In Sect. 2, we use perturbation methods to identify the dominant dynamics of the ST-PDE and P-PDE in parameter regimes that correspond to chemotaxis-dominated tip cell movement and small branching rates. We find that, in such regions of parameter space, both PDE models reduce to the same system at leading order within the domain interior. Numerical simulation confirms that pointwise differences between both solutions are relatively small under the conditions outlined above. In Sect. 3, we construct asymptotic solutions to the leading order systems, which are valid for early times. Numerical simulation demon-  strates that such analytic expressions are in good agreement with ST-PDE and P-PDE results. The approach taken in this paper may be used in the future to compare other nonlinear coarse-grained and/or phenomenological models for a given biophysical process.

Leading order model dynamics
The procedure we use to reduce the ST-PDE and P-PDE to their dominant, leading order dynamics is presented below. For clarity, the details are shown only for the ST-PDE: the analysis for the P-PDE follows mutatis mutandis, and is presented in Appendix A. We focus our analysis only on the domain interior: as we have already noted, the continuum models are simulated from t 0 = 0.2; at this time point, most of the tip cell mass has moved away from the left-hand boundary into the domain interior. Furthermore, both continuum models represent the migration of cells towards a TAF source, and not the subsequent infiltration of cells into the region producing TAF; this means that both models cease to be valid biological descriptions when the majority of the tip cell mass reaches the right-hand boundary at x = L. We therefore restrict our attention to solutions within the domain interior rather than near the boundaries. We make the following assumptions: (A1) All parameter values in the ST-PDE are equal to their counterparts in the P-PDE. (A2) The initial and boundary conditions for the ST-PDE and P-PDE are identical.
After substituting equation (7) into the ST-PDE stalk cell evolution equation, we simplify our analysis by reducing the number of parameters in the ST-PDE. We do this by recasting the dependent and independent variables with the following trans-formations: so that the time scale of interest is λ −1 and the length scale is √ D/λ. We write the TAF concentration as a function c(X , τ ) for purposes of clarity, as this allows us to avoid explicitly writing x in terms of the new independent variables. Equations (5)-(6) then transform to give with boundary and initial conditions given by In Eqs. (9)-(12), we have introduced the following dimensionless parameter groupings := √ Dλ χν , α := λ μ , β := a e μ a n λ = a e a n α .
We remark that, for this specific scenario, β is a constant because κ(x) and the TAF gradient are constant; in more general scenarios where the TAF gradient varies in space, however, β would be a function that also depends on the rescaled spatial variable. Equations (9)-(12) may be compared to the rescaled P-PDE (Eqs. (21)-(22) in Appendix A). Under the following assumption, we may further reduce the ST-PDE using perturbation methods (Bender and Orszag 1999;Hinch 1991;Verhulst 2005): (A3) Tip cell movement is dominated by chemotaxis and the branching rate is sufficiently small so that 0 < 1 and 0 Assumption (A3) is motivated by research that suggests chemotaxis up TAF gradients is the dominant mechanism by which tip cells migrate (Duran et al. 2017;Betz et al. 2016;Bowersox and Sorgente 1982;Carmeliet and Jain 2011). The relationship between and α follows from assumption (A1) and the relationship between the coarse-grained model parameters and those of its underlying discrete model (Simpson et al. 2009;Pillay et al. 2017). Namely, if k and h represent P-ABM parameters related to the response of tip cells to TAF and the spatial step size of the lattice, respectively, then the continuum parameters D and χ may be written as Substituting these expressions into the formula for , we find that = √ α/(2khν) = √ αΨ −1/2 . We will also assume for the remainder of this section that the ratio a e /α is bounded in the limit as α → 0, so that β ∼ O(1). In practice, this latter assumption may be relaxed: in Appendix B, we show that the leading order dynamics of the ST-PDE and P-PDE are identical for a wider range of values of β.
With 0 < 1, we seek asymptotic solutions for u(X , τ ) and w(X , τ ) of the form Substitution into Eq. (10), and equating coefficients of O( 0 ), leads to the following PDE, which describes the dynamics of w 0 (X , τ ): Ideally, the same procedure would be used to identify the leading order dynamics of the tip cell rate Eq. (9); however, this naive approach will fail because the equation contains a term of O( −1 ): in particular, we would find that the leading order tip cell solution gradient is 0, which clearly does not match our observations of the ST-PDE and P-PDE solutions in Fig. 1. We may circumvent this issue by making a change of variables: we let where X 0 is an arbitrary constant that ensures the right boundary location, y * = (λL/(χ ν)−τ − X 0 )/ , is strictly positive over the time interval of interest. Once again, we have written the steady state TAF concentration as a function C(y, τ ) for purposes of clarity. We remark that C(y, τ ) is still bounded between 0 and 1. Substitution of these transformations into Eq. (9) yields with boundary conditions  Parameter values: D = 10 −3 , χ = 0.4, a n = 1, a e = 0.0391, μ = 160 (this corresponds to = 10 −3/2 , α = 10 −3 , and β = 39.1). The PDEs were simulated on τ ∈ [0.2λ, 2λ], X ∈ [0, λ D ]. P-ABM solutions at τ = 0.2λ, column averaged in the y-direction, were used as initial conditions. For colours, we refer to the online article We seek a regular perturbation series solution for U of the form At leading order, we recover the following PDE describing the dynamics of U 0 : where the dynamics of W 0 are described by Eq. (13) (after a suitable change of variables). The boundary conditions for the leading order solution, found in the limit as which are homogeneous Dirichlet boundary conditions. We show in Appendix C that tip cell solutions to Eq. (16) with boundary conditions (17) are non-negative, provided the initial condition is also non-negative. We may therefore ignore the absolute value sign in Eq. (13). The resulting system is identical to that obtained from the P-PDE (see Appendix A). Figure 2 presents the numerical solution to the leading order system, with the tip cell solution transformed back into the independent variables X and τ , along with numerical results from the original ST-PDE. We observe good agreement between both sets of tip cell densities in Fig. 2a: at each time point shown, for instance, the largest pointwise difference between the two solutions is less than 2% of the maximum tip cell density. Similarly, we observe in Fig. 2b that the numerical solution to Eq. (13) is almost indistinguishable from that of the original ST-PDE when the leading order tip cell solution is used to calculate the stalk cell density rate of change. Indeed, at each time point shown the maximum pointwise difference between the two solutions is less than 2% of the maximum stalk cell density. These results support our earlier hypothesis that the leading order dynamics of the ST-PDE within the domain interior are accurately described by Eqs. (13) and (16) when 0 < 1 and 0 < α 1. Crucially, we show in Appendix A that, under assumptions (A1)-(A3) listed above, we may reduce the P-PDE to Eqs. (13) and (16) as well. This result holds regardless of the value of β, the non-dimensional parameter that depends on the rate of tip-to-sprout anastomosis (see Appendix B for details). In biological terms, this suggests that the ST-PDE and P-PDE are identical to leading order if chemotaxis dominates tip cell movement and if the branching rate is sufficiently small; additionally, the agreement between the two models does not depend on the rate of tip-to-sprout anastomosis. Since the higher-order terms in the tip and stalk cell asymptotic series will be small when 0 < 1 and 0 < α 1, this implies that solutions to the original ST-PDE and P-PDE models will be indistinguishable within the domain interior for such parameter regimes.
We have validated these conclusions with numerical results. Figure 3 presents numerical solutions of the original ST-PDE and P-PDE models listed in Table 1, using a parameter regime for which 0 < 1, 0 < α 1, and β ∼ O(1) (parameter values are listed in the figure caption; see Appendix E for details of the numerical methods). We observe in Fig. 3a that at each time point, the tip cell solution profile resembles a bell curve. We provide a possible explanation for this in Sect. 3, where we show that the leading order tip cell solution has a bell curve profile for early times, regardless of the value of β. We also see that the amplitude of the tip cell density increases gradually as the solution travels to the right. This most likely occurs because the tip cell branching rate increases near the right-hand boundary (the proliferation rate is proportional to the TAF concentration, which itself is linear with respect to X ).
The ST-PDE and P-PDE tip cell solutions appear to travel with similar speeds in Fig. 3a. In fact, further investigation shows that the P-PDE speed is within 0.5% of the ST-PDE wave speed, and that in both models the waves accelerate gradually over time (see Appendix E for details on the numerical methods used to arrive at this conclusion).
The ST-PDE tip cell solutions appear to be indistinguishable from those of the P-PDE. We investigated this further by computing the largest pointwise difference between the ST-PDE and P-PDE tip cell solutions over the spatial domain and normalising that value by the maximum tip cell density at each time point (results are presented in Fig. 3c). We observe that the maximum relative pointwise difference between the two sets of solutions lies within the approximate range [0.015, 0.035] (or 1.5%-3.5% of the maximum tip cell density) for the time points shown. Figure 3b shows the corresponding ST-PDE and P-PDE stalk cell densities for this parameter regime. We see that the P-PDE solution is greater than that of the ST-PDE, even though the solutions to both models appear to have similar profiles. Nevertheless, , a e = 0.0391, a n = 1, μ = 160 (this corresponds to = 10 −3/2 , α = 10 −3 , β = 39.1). The PDEs were simulated on The P-ABM solutions at τ = 0.2λ, column averaged in the y-direction, were used as initial conditions for the PDE models. For colours, we refer to the online article the size of differences between the ST-PDE and P-PDE solution is small: we calculated the maximum relative pointwise difference between stalk cell solutions using a similar method to the one described above for tip cell solutions (results not shown). In this case, the relative difference between the two sets of stalk cell solutions lies within the approximate range [0.01, 0.02] (1%-2% of the maximum stalk cell density) for the time interval considered here.
In Fig. 3a, the tip cell density initially decreases at a relatively rapid rate. We investigated this behaviour further by creating a log-log plot of the maximum tip cell density against time τ in Fig. 3d. The graph resembles a straight line until τ ≈ 0.4, after which point it begins to curve upwards. In fact, we found that a line with slope −1.02 approximates the graph with less than 2% error until τ ≈ 0.4 (the equation for the line was computed using least squares regression). We conclude that, for this parameter regime, the maximum tip cell density is approximately proportional to τ −1 at early times. We explain this observation in Sect. 3, where we find that for this parameter regime there exist self-similar solutions to the leading order ST-PDE and P-PDE dynamics that are initially proportional to τ −1 .
We confirmed that increasing the value of in the ST-PDE and P-PDE increases the magnitude of differences between their solutions. We numerically simulated the PDE models for cases in which the random motility coefficient, D, was increased from 10 −3 to 10 −1 , and in which the chemotactic sensitivity of cells, χ , was decreased from 0.4 to 0.04; in both cases increases from about 0.03 to approximately 0.3 but α and β are the same as in Fig. 3 (see Figures 1-2 of the Supplementary Material). In both situations, we were able to visually distinguish between the tip and stalk solutions of the two models. Further numerical calculation also revealed that, for such parameter regimes, the relative differences between tip and stalk cell solutions of the two models were larger than the differences shown in Fig. 3.
We have also investigated the case in which the values of α and increase simultaneously by changing the value of λ, the parameter in the ST-PDE and P-PDE related to the branching rate, from 10 −3 to 0.0625 (see Figure 3 of the Supplementary Material). Interestingly, for this parameter regime the differences between the ST-PDE and P-PDE tip/stalk solutions grow noticeably large only over long time periods: the two sets of solutions are indistinguishable from each other at early time points. This apparent discrepancy from our analysis likely arises because the branching rate is proportional to λ and to the TAF concentration. Near the left-hand boundary, the TAF concentration is small because it is linear with respect to X . In such regions, the overall branching rate remains low even when the value of λ increases, and it follows that solutions travelling in such regions are less likely to be affected by changes in λ. As the tip cell density travels to the right, however, the TAF concentration and branching rate both increase and become more important to solution behaviour. Thus if only the value of λ were to increase, then differences between ST-PDE and P-PDE solutions would grow only over long time periods, when the tip cells have invaded regions with sufficiently large TAF concentrations.
Our analysis suggests that β does not play an important role at leading order, hence its value is not expected to affect the degree of similarity between the ST-PDE and P-PDE solutions. We confirmed this conclusion via numerical simulation, by increasing the value of a e in the ST-PDE and P-PDE from 0.0391 to its maximum value of 1 (see Figure 4 of the Supplementary Material; this corresponds to increasing β from 39.1 to 10 3 but keeping and α the same values as in Fig. 3).
We have thus found that numerical simulation corroborates the conclusions from our analysis in Sect. 2. Namely, we have found that increasing the value of = √ Dλ/(χ ν) and/or α = λ/μ increases differences between the ST-PDE and P-PDE solutions. However, changing the value of β does not affect whether the two sets of numerical solutions appear indistinguishable from each other. Hence if 0 < 1 and 0 < α 1, then solutions for the ST-PDE are in good agreement with those of the P-PDE. This corresponds to biological scenarios in which chemotaxis dominates tip cell movement and sprout branching is rare.

Self-similar leading order tip cell solutions can exist for small time periods
In this section, we construct asymptotic solutions to the leading order tip cell dynamics derived in Sect. 2. Specifically, we show that Eq. (16), which describes the leading order dynamics of the ST-PDE and P-PDE within the domain interior, admits selfsimilar solutions for early times (i.e., in the limit as τ → 0). The results in Fig. 3d motivate our search for such solutions since we observe that, for early times, the maximum values of the tip cell solutions to the original ST-PDE and P-PDE models are proportional to the same fixed power of τ . We present the derivation of such self-similar solutions in Appendix D and list the main results of that analysis here. In the case where β ∼ O(1), for instance, we find that the leading order tip cell solution U 0 (y, τ ) may be written in terms of the self-similar expression τ −1 U 0 (z), where z = τ −1/2 (y − σ (τ )τ ) and σ (τ ) is the wave speed at time τ . Using the principles of dominant balance, we find that in the limit as τ → 0 the ordinary differential equation (ODE) describing the dynamics of U 0 (z) is which describes the self-similar tip cell solution profile. The terms that appear in Eq. (18) represent tip-to-tip anastomosis and the movement of tip cells. Terms due to branching and tip-to-sprout anastomosis are proportional to τ (see Appendix D for details, these correspond to the terms CU 0 and U 0 W 0 in Eq. (16)), and hence are ignored in the limit as τ → 0. Since these terms will become non-negligible for sufficiently large values of τ , we also conclude that this self-similar solution breaks down as τ increases because it becomes impossible to reduce the dynamics for U 0 to one independent variable. We thus reason that this self-similar solution exists when terms representing branching and tip-to-sprout anastomosis are negligible compared to those describing tip-to-tip anastomosis. To our knowledge, there is no closed-form solution to the boundary value problem given by Eq. (18) that is non-trivial (i.e., not U 0 (z) ≡ 0). However, a non-trivial numerical solution to the boundary value problem can be computed; it is presented in Fig. 4a. We see that the solution resembles a bell curve, which is consistent with our earlier observations of numerical tip cell solutions. We further compared the similarity solution to results for the ST-PDE and P-PDE by transforming the function in Fig. 4a from an expression in terms of the independent variable z into one in terms of the unscaled variable X : if we redefine the wave speed as φ(τ ) := 1/ + σ (τ ), then this transformation is given by  which follows from the relationships listed above and in Sect. 2 (see Appendix E for details on the numerical methods used to estimate φ(τ ) and X 0 ). The transformed selfsimilar solution is plotted in Fig. 4b over several values of τ for a parameter regime in which = 10 −3/2 , α = 10 −3 , β = 39.1, and λ = 0.16; the numerical solutions to the original ST-PDE and P-PDE models are also plotted for comparison. We observe in Fig. 4b that the self-similar solution is in good agreement with the numerical solutions for the ST-PDE and P-PDE for small values of τ (the largest pointwise difference between the self-similar solution and the numerical results is on the order of 10% of the maximum tip cell density). As τ increases, however, the self-similar solution begins to underestimate the ST-PDE and P-PDE solutions. This observation is consistent with our earlier argument that self-similar solutions break down as τ increases. The results shown in Fig. 4 thus provide evidence that the leading order system admits a similarity solution when τ is sufficiently small, and that its profile is well described by Eq. (18).
At this point, we have only presented self-similar solutions to the leading order inner tip cell density in the case where β ∼ O(1). However, similarity solutions also exist for a wider range of values of β (this can be shown using a similar procedure to the one presented in Appendix D). When β 1, for instance, the self-similar solution is identical to the one described above. This makes sense, since we show in Appendix B that terms representing tip-to-sprout anastomosis in the PDEs describing the leading order dynamics are negligible for this parameter regime; thus, this result is consistent with our conclusion that Eq. (18) is an appropriate description of the leading order solution when terms describing branching and tip-to-sprout anastomosis are negligible compared to those for tip-to-tip anastomosis.
When β 1, however, the self-similar solution is different to the one described by Eq. (18). This is because we find in Appendix B that, for this parameter regime, terms corresponding to tip-to-tip anastomosis are negligible in the leading order tip cell solution. A similar procedure to the one listed in Appendix D demonstrates that in the limit as τ → 0, the ODE describing the self-similar solution profile is given by where z = τ −1/2 (y − σ (τ )τ ) is the same similarity variable as before and A is a constant such that U 0 (y, τ ) = τ A U 0 (z). The terms in this equation represent tip cell movement only: terms describing tip-to-sprout anastomosis and branching vanish as τ → 0. We conclude that the self-similar solution described by Eq. (19) is applicable in cases where terms that represent branching, tip-to-tip anastomosis, and tip-to-sprout anastomosis are negligible compared to those for tip cell movement. The value of A is uniquely determined by solving the ODE, satisfying the two boundary conditions, and enforcing the constraint that the solution must be non-negative. This sets A = −1/2, and the solution to Eq. (19) is where C 1 is a non-negative constant. The solution is a Gaussian function and thus has a bell curve profile. We have verified that it can be in good agreement with the ST-PDE and P-PDE numerical solutions for early times and parameter regimes for which 0 < 1, 0 < α 1, and β 1 (see Figure 5 of the Supplementary Material).

Discussion and conclusion
In this article, we have analysed two continuum models for angiogenesis: the snail-trail and Pillay PDE systems. The P-PDE is a nonlinear mean-field model that describes the coarse-grained behaviour of a discrete, rule-based ABM (Pillay et al. 2017). By contrast, the snail-trail model is a simpler system constructed from phenomenological arguments (Byrne and Chaplain 1995;Balding and McElwain 1985;Orme and Chaplain 1997;Flegg et al. 2020;Connor et al. 2015). Despite such differences, asymptotic solutions to the two models are identical at leading order within the domain interior under parameter regimes in which chemotaxis dominates tip cell movement and sprout branching is rare. Hence in such biological situations we anticipate that solutions to the ST-PDE will be in good agreement with P-PDE results; we confirmed this with numerical simulation.
In Sect. 3 we demonstrated that the system describing the leading order dynamics of the ST-PDE and P-PDE admits self-similar solutions for early times. Differences in numerical results between the similarity solution, ST-PDE, and P-PDE are relatively small, which indicates that higher order terms in the asymptotic expansions are indeed negligible under parameter regimes in which tip cell movement is dominated by chemotaxis and branching is rare. Our analysis of leading order solutions also provides insight into some observations of the numerical results. For instance, the leading order similarity solution has a bell curve shape regardless of the value of β, the nondimensional parameter that is proportional to the rate of tip-to-sprout anastomosis; this may explain why tip cell solutions evolve to such a profile over time.
There are several possible extensions to our work in this article. We limited our search of asymptotic solutions to those only within the domain interior, as we reasoned that they were most relevant to the numerical results given the models and initial conditions. However, this does not exclude the possibility that boundary layers may exist: indeed, numerical results (not shown) reveal that stalk cell solutions of the ST-PDE and P-PDE are different to each other near the left-hand boundary. Thus, one extension of our analysis is to identify potential boundary layers in the leading order ST-PDE and P-PDE solutions, and match them to the dynamics derived in Sect. 2. We also made a quasi-steady state approximation for the TAF dynamics, which simplified our analysis by making the TAF field linear. TAF dynamics may be more complicated in reality, however, and can play a more significant role in the migratory dynamics: cell-induced gradients, for instance, can provide time-dependent guidance cues for tip cells, and will generate spatial heterogeneities in the TAF gradient that were not considered here (Tweedy et al. 2020).
The approach taken in this paper may be generalised to study the relationship between the ST-PDE, P-PDE, and other nonlinear continuum models for angiogenesis. As a particular example, we cite a coarse-grained system from Bonilla et al. (2016), which has been derived from a discrete model that relies on stochastic differential equations to update cell and vessel locations (Bonilla et al. 2014). The continuum model of Bonilla et al. is different from the ST-PDE and P-PDE, but its solutions nevertheless have similar behaviours to the ones investigated in this article: for instance, the leading order tip cell solution in the Bonilla et al. PDE model has been shown to have a profile that resembles a bell curve, and is a travelling wave (Bonilla et al. 2016(Bonilla et al. , 2020. It would be interesting in future work to apply our approach to these other continuum angiogenesis models, as it would provide analytic insight into when solutions from different angiogenesis models can be distinguished from each other, and when simple phenomenological models will be accurate representations of more complicated systems. Although we identified parameter regimes for which two continuum models are expected to be in good agreement with each other, one open question we did not address in this work is when the continuum models will also closely approximate solutions of discrete models for angiogenesis (which are closer representations of the true biology). Some insight into this question may be gained by considering the derivation of the P-PDE, since this continuum model describes the mean-field behaviour of an ABM and has the same leading order solution as the ST-PDE. In particular, mean-field models such as the P-PDE are appropriate representations of discrete models when spatial correlations between agents can be neglected. In the case of on-lattice ABMs, for instance, this means that the average occupancy of a given lattice site must be independent of the occupancy of its neighbours (Simpson and Baker 2011;Pillay et al. 2017Pillay et al. , 2018. However, this mean-field assumption will be violated in most discrete angiogenesis models because of volume exclusion rules due to anastomosis, and the effect of this violation will become more noticeable when many cells are located within a given spatial region (Simpson et al. 2010). We thus expect the P-PDE (and, by extension, the ST-PDE) to be a poor description of discrete angiogenesis results when tip cells branch often. Figure 5 presents column averaged P-ABM solutions for a parameter regime in which the branching probability P p has been increased from 10 −3 to 5 × 10 −2 , along with ST-PDE and P-PDE solutions (several continuum parameters have been fitted to the ABM results with a nonlinear least squares method, see the figure caption for their values). We observe that the two continuum models are indistinguishable and hence are in good agreement (the largest pointwise difference between the ST-PDE and P-PDE tip cell solutions at each time point, for example, is less than 4% of the maximum tip cell density). However, the discrepancies between the continuum models and the P-ABM distributions are larger, because we can visually distinguish between the continuum and discrete results. Such results confirm that the continuum models will be increasingly poor descriptions of discrete ABM results when the branching rate increases, even though the two continuum models may still be in good agreement with each other.
Surprisingly, further inspection of Fig. 5 shows that the continuum models do not completely fail at capturing some aspects of the discrete solutions. For example, the shape of the continuum solutions resembles those of the column averaged discrete results, despite the fact that the sizes of their amplitudes are different. Additionally, the continuum models qualitatively capture the experimentally observed brush-border effect, since the tip and stalk cell densities increase as they move closer to the righthand boundary at x = 1. Most importantly, we observe that the leading edge of tip and stalk cell solutions travel with almost identical speeds (they are within 3.5% of each other). We conclude that the continuum models are able to capture some aspects of the underlying discrete solutions, such as the speed of invasion and their qualitative behaviour, even though the calculated densities do not entirely match with the discrete results.
This raises the question of how one should quantify the degree of agreement between discrete and continuum results. For example, "good agreement" could be characterised by sufficiently small absolute pointwise differences between the discrete and continuum solutions, or by small differences between their wave speeds. Clarifying what it means for discrete and continuum models to be in good agreement may provide insight into new metrics that could indicate both the extent of angiogenic progression and when it would be appropriate to apply simple continuum models to study and analyse in silico (and eventually in vivo) angiogenesis data.
We anticipate that some of our results from this article also hold in a discrete setting due to the relationship between the P-PDE and its underlying ABM. In particular, we anticipate that for parameter regimes in which tip cell branching is rare and cell movement is dominated by chemotaxis, solutions to the two continuum models will   N (x, t), and b stalk cell density E(x, t), given by the ST-PDE, P-PDE, and P-ABM at times t = 0.2, 0.4, . . . , 2. One thousand realisations of the P-ABM, with branching parameter P p = 5 × 10 −2 , were column averaged in the y-direction to generate the results (all other ABM parameter values are the same as those listed in Appendix E). Column averaged P-ABM solutions at t = 0.2 were used to initialise the two PDE models, which were simulated on the interval . For colours, we refer to the online article be in good agreement with ensemble averages of the P-ABM. If this is true, then this suggests that several summary statistics, such as the number of branches and the average tip cell displacement, could be used to anticipate when discrete solutions would be captured well by continuum models. It would be interesting in future work to evaluate how well such statistics predict the level of similarity between P-ABM and ST-PDE solutions, as this will inform us of when simple phenomenological systems will be good descriptions of more complex discrete angiogenesis data, potentially including those obtained from in vitro and in vivo approaches.
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/.

A Derivation of leading order P-PDE dynamics
Here we derive the leading order dynamics of the P-PDE, using the same procedure and assumptions listed in Sect. 2.
We recast the independent and dependent variables of the P-PDE in terms of the transformations given by Eq. (8), so that N (x, t), E(x, t), and C(x, t) map to the functions u(X , τ ), w(X , τ ), and c(X , τ ), respectively. This leads to the equations which are subject to the initial and boundary conditions given by Eqs. (11)-(12) in Sect. 2. Under assumption (A3) of Sect. 2, we let 0 < 1, 0 < α = 2 Ψ 1, where Ψ is a constant of O(1). We may thus further reduce the P-PDE using perturbation methods. We will assume for the remainder of this appendix section that β ∼ O(1) (we show in Appendix B that the leading order dynamics are the same for the ST-PDE and P-PDE, even for different values of β).
Using the same procedure as outlined in Sect. 2, we find that the dynamics of the leading order stalk cell solution, w 0 (X , τ ), are given by The dynamics of the P-PDE leading order tip cell solution are found by transforming the independent and dependent variables in Eq. (21) as where X 0 is an arbitrary constant, so that u(X , τ ), w(X , τ ), and c(X , τ ) map to functions U (y, τ ), W (y, τ ), and C(y, τ ). This leads to the equation We substitute the perturbation series U (y, τ ) ∼ U 0 (y, τ ) + U 1 (y, τ ), and identify the leading order dynamics by equating terms of O( 0 ). This leads to the leading order P-PDE tip cell equation where the dynamics describing the leading order stalk cell solution, W 0 (y, τ ), are given by Eq. (23) (after a suitable transformation of variables). Crucially, the resulting inner system is equivalent to that of the ST-PDE in Sect. 2, since we show in Appendix C that the leading order tip cell solutions to the ST-PDE are non-negative. We have confirmed that numerical solutions to the above inner system are in good agreement with results from the full P-PDE (see Figure 6 of the Supplementary Material).

B Leading order dynamics for different values ofǏ
n Sect. 2 and Appendix A, we assumed β ∼ O(1) when proving that the ST-PDE and P-PDE are identical at leading order. We may, however, relax this condition: this section presents the procedure for reducing the ST-PDE to leading order for the cases in which β 1 and β 1 (the analysis for the P-PDE follows naturally, and is omitted for clarity).
When β 1, we may assume that β = β −m , where m > 1 and β is a constant of O(1). To derive the leading order dynamics, we rescale the dependent variable u so that u = mū and u(X , τ ) →ū(X , τ ). Substitution into Eqs. (9)-(10) yields From this point, the procedure for deriving the outer and inner systems is the same as in Sect. 2. We find that the dynamics governing the leading order stalk cell density, w 0 (X , τ ), are given by Similarly, the leading order tip cell dynamics,Ū 0 (y, τ ), may be identified using the same transformations presented in Sect. 2. They are given by In all cases, u(X , τ) and U (y, τ) represent the tip cell density, w(X , τ) and W (y, τ) the stalk cell density, and c(X , τ) and C(y, τ) the TAF concentration. We refer to Sect. 2 and Appendices A and B for the derivation of these systems where the dynamics of W 0 (y, τ ) are given (after a suitable transformation of variables) by Eq. (28). The boundary conditions are given at leading order by lim y→±∞Ū 0 (y, τ ) = 0.
A similar proof to the one in Appendix C shows that the leading order tip cell solution is non-negative, provided its initial condition is non-negative. Therefore we may drop the absolute value sign in Eq. (28). Crucially, the resulting leading order system is identical to the one derived for the P-PDE. For the case in which β 1, we rewrite the parameter in terms of a power of , such that β = β m , where m > 0 and β is a constant of O(1). We rescale the dependent variable w in the ST-PDE, such that w = mw and w(X , τ ) →w(X , τ ). Substitution into Eqs. (9)-(10) leads to the system We deduce that the dynamics of the leading order stalk cell solution,w 0 (X , τ ), are given by while the dynamics of the leading order tip cell solution are given by where y = X − τ +X 0 , subject to the homogeneous Dirichlet boundary conditions If the initial condition for the leading order tip cell solutions is non-negative, then a similar proof to the one outlined in Appendix C shows that their solutions are nonnegative at all other times. Thus we may neglect the absolute value signs in the leading order stalk cell evolution equation. The resulting inner and outer systems are identical to those given by the P-PDE. Table 3 summarises all of the possible outer and inner systems for the ST-PDE and P-PDE when 0 < 1 and 0 < α 1.

C Leading order ST-PDE solutions are non-negative
We show in this appendix that the solution to the leading order ST-PDE tip cell dynamics given by Eq. (16) are non-negative, provided the initial condition is also non-negative. We will assume that solutions are at least twice differentiable. For sake of contradiction, assume that there exists a point y * ∈ R such that at time τ 1 > 0, U 0 (y * , τ 1 ) < 0. Since the PDE solution is twice differentiable, it is continuous, so there must have existed some time 0 ≤ τ * < τ 1 such that U 0 (y * , τ * ) = 0 and ∂U 0 ∂τ (y * , τ * ) < 0. We also conclude that U 0 (y, τ * ) ≥ 0 for all other values of y, since y * is the first point at which the solution becomes negative.

D Derivation of self-similar tip cell solutions
In this appendix, we show that the leading order tip cell solution is self-similar in the limit as τ → 0. We focus on the case for which β ∼ O(1); a similar procedure can be used to derive self-similar solutions for a wider range of values for β, using the inner leading order dynamics derived in Appendix B.
In addition to the assumptions made in Sect. 2, we will simplify our analysis by assuming that the following are true: (B1) The leading order tip cell solution dynamics are given by Eq. (16) with boundary conditions (17). (B2) The TAF concentration C(y, τ ) is bounded between 0 and 1, as is the case in Sect. 2. (B3) The initial leading order tip and stalk cell densities are non-negative and bounded. (B4) The wave speed of leading order tip cell solutions, σ (τ ), may depend on time but is bounded such that σ (τ )τ and σ (τ )τ → 0 as τ → 0. This assumption is inspired by our numerical results, which suggest that solutions to the ST-PDE and P-PDE travel short distances in small time periods but have wave speeds that vary over time (see Fig. 6). We seek self-similar solutions by substituting into Eq. (16) the expression where A and B are constants to be determined, σ (τ ) is the tip cell wave speed at time τ , and z = τ B (y − σ (τ )τ ) is the self-similar variable. We recast the equation in terms of the independent variables z and τ . To avoid explicitly writing the TAF concentration as a function of z, we write C(y, τ ) as c (z, τ ). Similarly, we use W 0 to denote the leading order stalk cell density in terms of z and τ . These substitutions lead to the equation The values of A and B are chosen to eliminate as many powers of τ in Eq. (36) as possible, while simultaneously ensuring that no term blows up to infinity as τ → 0. This uniquely sets A = −1 and B = −1/2. However, the resulting equation does not describe self-similar solutions in general. We conclude that self-similar solutions do not exist for sufficiently large values of τ . However, in the limit as τ → 0 one can recover an ODE for U as a function of z. This is because the TAF concentration, wave speed, and stalk cell densities are all bounded: since such quantities are multiplied by positive powers of τ in Eq. (36), they vanish when τ → 0. Such simplifications lead to the ODE (18).

E Numerical methods
Initial conditions for the ST-PDE and P-PDE were constructed from ensemble averages of the P-ABM, which were generated using the following algorithm from Pillay et al. (2017). ABM parameter values are listed in Table 4.
Choose Δt, the time step, and t final , the time at which to terminate the ABM solution.
where c(x i , y j ) is the TAF concentration at location (x i , y j ). (iii) If the tip cell moves into a site occupied by one (or more) tip cells, then tip-to-tip anastomosis occurs and the tip cells are removed. (iv) Otherwise, if tip cell moves into a site occupied by a stalk cell that does not belong to the same sprout, then tip-to-sprout anastomosis occurs and the tip cell is removed.

End While Loop
One thousand realisations of the P-ABM were column-averaged in the y-direction at t = 0.2 to set up initial conditions for the PDE models. We initialised the PDEs at time t = 0.2, rather than t = 0, as previous investigations have shown that the P-PDE can become ill-posed if it is initialised at earlier times (Pillay et al. 2017(Pillay et al. , 2018. The 1D PDE models were solved using the method of lines (Schiesser 1991). Using grid points with spatial step size Δx = 1/200, spatial derivatives were discretised using central finite difference formulae. The systems were solved with MATLAB's ode15s solver, which implements adaptive time-stepping for stiff ODE problems (Shampine and Reichelt 1997).
In Sects. 2 and 3, we transformed solutions of the ST-PDE and P-PDE from functions of the form N (x, t) and E(x, t) to u(X , τ ) and w(X , τ ), respectively, so as to facilitate comparison with our analysis. This is done using the formulae listed in Eq. (8).
The solution to Eq. (18), presented in Fig. 4a, was computed numerically using the MATLAB library Chebfun, which applies spectral methods and Newton's method to solve nonlinear boundary value problems (Driscoll et al. 2014). To approximate the infinite domain, the solution was computed on the interval z ∈ [−100, 100], which is two orders of magnitude larger than the region over which the solution gradient varies the most. The numerical solution was computed using the initial guess U 0 (z) = exp(−z 2 ).
The quantities X 0 and φ(τ ), which are used to transform the self-similar solution into the independent variable X in Fig. 4b, are calculated in the following manner: the wave speed φ(τ ) is taken to be the derivative of X (τ ), a function that describes the tip cell wave front location at time τ . We define the wave front location as the mean location of the tip cell solutions: In the case where X (τ ) is approximately linear, φ(τ ) is estimated as the slope of the line that minimises the squared error with the data of X (τ ) versus τ , and the value of X 0 is the y-intercept of this line of best fit. Figure 6 shows a plot of X (τ ) against τ for the same parameter regime and initial conditions listed in Fig. 3. We observe that the ST-PDE and P-PDE models have roughly identical wave speeds for the time points considered here, and that the ST-PDE and P-PDE tip cell solutions gradually accelerate over time. For early times, however, the two solutions appear to travel with constant speed: we confirmed this by plotting a straight line that minimises the squared error with the numerical solutions (the line was fit to data in the range τ ∈ [0.2λ, 5λ], where λ = 0.16). We see in Fig. 6 that the line of best fit is a good approximation to the numerical results for sufficiently small values of τ (for τ < 0.8, the R 2 -value is greater than 0.99). We conclude that for τ ≈ 0, the wave speed φ(τ ) is approximately constant; its value is given in Table 2 of Sect. 3.

F Derivation of snail-trail model corrective factor
In this appendix, we provide a brief overview of our derivation of κ(x), the corrective factor in the snail-trail model given by Eqs. (5)-(6). For more details on the derivation, we refer to the original article by Martinson et al. (2020).
We let κ be a non-dimensional factor in the snail-trail model that corrects for neglected vessel production in directions other than that of the migrating front. Since a new vessel is itself made up of stalk cells, κ is related to κ through the relationship κ = κ/h EC , where h EC represents the length of a typical stalk cell. Either factor must be included in the snail-trail modelling framework because the tip cell flux, a net quantity, tends to underestimate the total vessel/stalk cell density rate of change. Thus, if J net measures the average net tip cell flux in the x-direction within a time step Δt and spatial interval h, then κ|J net | is equal to the (true) density of new vessels that are created on average within this interval (note that we have absorbed any dependence of h and Δt into κ). In order to derive a function for κ (and, by extension, κ), we examine how the above statement translates into a discrete setting.
We consider the P-ABM framework as a discrete representation of tip and stalk cell movement following a snail-trail assumption. In the P-ABM, tip cell movement is a 2D biased random walk towards increasing concentrations of TAF. The ABM is assumed to have N lattice sites that are equally spaced with step size h, where h represents the length of a typical cell. We assume, for sake of simplicity, that a tip cell length is equivalent to that of a stalk cell (so h EC = h). Just as for the continuous snailtrail model, we assume that the TAF concentration C(x) is prescribed and independent of time t and the spatial variable y.
Within a time step Δt, tip cells are chosen to move with a constant probability P m . Once chosen to move, the direction in which a tip cell travels is selected according to the following probabilities: if we denote the probability of moving left as P x − and right as P x + , then these values are defined at x i ∈ (0, 1) as where for 0 < i < (N − 1), For example, if g x = −1/2, then P x − = 3/8 and P x + = 1/8 (i.e., the tip cell is more likely to move in the left direction). Similar formulae are used to determine the movement of tip cells in the y-direction. The parameter value k in Eq. (39) ensures that no numerator in Eq. (38) can become negative (so that |g x (x)| ≤ 1); the value of k does not vary with respect to location. In the P-ABM, a discrete version of the snail-trail assumption is incorporated by having stalk cells proliferate in the space left empty by a moving tip cell. In other words, any new vessel production within a time step Δt is equal to the number of times tip cells move, since stalk cells will proliferate to occupy the resulting empty space.
This leads us to immediately identify a discrete analogue of |J net |, the magnitude of the tip cell net flux in the x-direction: it is simply the expected net number of jumps that tip cells make from a lattice site in the x-direction (we use expected values here because of the stochastic nature of the discrete model). We may then interpret the value of κ as the total amount of new vessel density produced for every unit of net tip cell movement in the x-direction.
We now use the rules of the P-ABM to compute the expected number of jumps that tip cells make in the x-direction. We define X R and X L , to be random variables that measure the total number of rightward and leftward jumps that originate from lattice point x i , respectively. Since cells may only move in these two directions, we model these random variables using multinomial distributions. We also define the random variable X net = X R − X L , which measures the net number of rightward jumps made from lattice site (x i , y j ). By the argument above, the net tip cell flux in the x-direction is defined by J net := E[X net ]. Its magnitude is given as |J net |.
Since we have reasoned that κ|J net | is the total vessel density produced (or, equivalently, total jumps that occur in the x-direction) from a given lattice site within the time step Δt and spatial interval h, it follows that when the value of |J net | is normalised to 1, κ becomes equal to the total number of jumps in the x-direction from lattice site (x i , y j ). We exploit this information to calculate the expected values of our random variables: for example, the probability of executing m rightward jumps becomes where κ m is the binomial coefficient. A similar equation holds for X L . The expected values of X R and X L are thus given by when |J net | is normalised to 1, so that Since Eqs. (40) From this equation, it is clear that κ is non-negative.
We may apply Taylor's theorem to simplify Eq. (42): assuming that the average cell length is sufficiently small so that 0 < h 1, we may write g x as Substitution of this expression into Eq. (42) yields hence κ(x) is inversely proportional to the magnitude of the local TAF gradient. It is possible to further transform Eq. (43) so that κ(x) is in terms of the continuum parameters D and χ . This simplification follows from the relationship between biased random walk models and advection-diffusion equations (Simpson et al. 2009;Pillay et al. 2017). Namely, if P m denotes the probability that a tip cell moves within a time step Δt, then the relationship between the discrete and continuum parameters is under the assumption that the above limits exist, are non-zero, and are finite. Substituting Eq. (44) into Eq. (43) yields where we have defined μ := P m /Δt. From here, it is straightforward to derive Eq. (7) by using the relationship κ(x) = κ(x)/h EC = κ(x)/h.