{'fraction_non_alphanumeric': 0.03931700966401282, 'fraction_numerical': 0.01923789494767413, 'mean_word_length': 4.693044469783352, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 12, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A novel nonlocal partial differential equation model of endothelial progenitor cell cluster formation during the early stages of vasculogenesis

Neovascularisation is essential for tissue development and regeneration, in addition to playing a key role in pathological settings such as ischemia and tumour development. Experimental findings in the past two decades have led to the identification of a new mechanism of neovascularisation, cluster-based vasculogenesis, during which endothelial progenitor cells (EPCs) mobilised from the bone marrow are capable of bridging distant vascular beds in a variety of hypoxic settings in vivo. This process is characterised by the formation of EPC clusters during its early stages and, while much progress has been made in identifying various mechanisms underlying cluster formation, we are still far from a comprehensive description of such spatio-temporal dynamics. In order to achieve this, we propose a novel mathematical model of the early stages of cluster-based vasculogenesis, comprising of a system of nonlocal partial differential equations including key mechanisms such as endogenous chemotaxis, matrix degradation, cell proliferation and cell-to-cell adhesion. We conduct a linear stability analysis on the system, solve the equations numerically, conduct a parametric analysis of the numerical solutions of the 1D problem to investigate the role of underlying dynamics on the speed of cluster formation and the size of clusters, and verify the key results of the parametric analysis with simulations of the 2D problem. Our results, which qualitatively compare with data from in vitro experiments, elucidate the complementary role played by endogenous chemotaxis and matrix degradation in the formation of clusters, and they indicate that previous approaches to the nonlocal modelling of cell-to-cell adhesion, while they capture the aggregating effect of cell-to-cell adhesion, are not sufficient to capture its stabilising effect on clusters, and new continuum cell-adhesion modelling strategies are required.


Introduction
The formation of new vascular networks is essential for tissue development and regeneration, as a functional vasculature is critical for tissue homeostasis. It is responsible for the delivery of oxygen and nutrients as well as the disposal of waste products. In addition, neovascularisation of local tissue is critical in a variety of pathological processes, among which are retinopathy, tumour growth, wound healing and soft-tissue ischemia. The development of a vascular network in localised solid tumours is particularly well known to promote further tumour growth and metastases. Therefore a better understanding of the mechanisms governing neovascularisation can help improve current therapeutic strategies, as well as identify new ones.
The two fundamental processes of neovascularisation are known as vasculogenesis, the de novo formation of a new vascular network, and angiogenesis, the formation of new blood vessels from pre-existing ones. While the distinction between these two processes has been clear since the early days of the field of study of vascular development [91], for a long time the term vasculogenesis was only associated with endothelial progenitor cells (EPCs)-mediated embryonic vasculogenesis, while blood vessels formation in adult organisms was believed to be formed predominantly, if not exclusively, via angiogenesis [88]. In the past 20 years great progress was made in understanding neovascularisation processes, in particular with the discovery of bone marrow-derived EPCs and postnatal vasculogenesis [9,11,10,56,83,88,101], for which a specific mechanism, here referred to as clusterbased vasculogenesis, has recently been proposed [16]. This newly identified mechanism of neovascularisation has been observed only in hypoxic settings in vivo, from development [84] to soft-tissue ischemia [109] to tumour vascularisation [115], and has been studied in detail by Blatchley et al. [16] in hypoxic gradients in vitro.

Cluster-based vasculogenesis
During this process, EPCs residing in the bone marrow are recruited in the circulation to the hypoxic site [10,101], where they first organise into clusters. They then proceed to sprout from these clusters and, in combination with angiogenesis, form the new vascular network merging with pre-existing vasculature. The term "clusterbased" vasculogenesis was therefore proposed by Blatchley and coworkers [16] to highlight the key step of the process which well distinguishes it from classical "single-cell" vasculogenesis, during which sparse mature endothelial cells (ECs) directly reorganise into a vascular network [98].
While EPCs recruitment to angiogenic sites has been well studied, we are still far from understanding how clusters form and what drives sprouting from these clusters. Akita et al. [1] reported that EPC differentiation, secretion of angiogenic factors -such as Vascular Endothelial Growth Factor (VEGF) -and migration were enhanced by hypoxic conditioning, resulting in increased EPC cluster formation in hypoxia in vitro. Nevertheless, it was not until recently that Blatchley and coworkers proposed, together with the name here used, a two-step mechanism of cluster formation and stabilisation in hypoxic environments [16], illustrated in Figure 2. The proposed mechanism relies on the observation that in highly hypoxic environments the production of matrix metalloproteases (MMPs) is upregulated by the presence of reactive oxygen species ( Figure 2B). This results in significant extracellular matrix (ECM) degradation responsible for EPCs aggregation into clusters, which are then stabilised by cell-to-cell adhesion ( Figure 2C). EPC clusters were observed between 24 and 48 hours from the beginning of the experiments, presenting diameters in the range 100 − 400 µm-see Figure 1. Finally, after 48 hours, cell-to-matrix interactions lead to sprouting and network formation, particularly fostered by highly viscoelastic -where not degraded -ECM ( Figure 2D). The final network, observed after 72 hours, was reported to be up to 500 µm in length, larger than previously observed ones formed via angiogenesis or single-cell vasculogenesis. Overall, cluster-based vasculogenesis has been shown to be of great importance in highly hypoxic tissues, as EPC clusters form in these regions far from pre-existing blood vessels and, after sprouting, function as a bridge for distant vascular beds, allowing for neovascularisation in physiological and pathological settings in which angiogenesis alone would not have sufficed.
While this is a great step forward in the understanding of cluster-based vasculogenesis, further work is required to reach an exhaustive comprehension of the process, as well as to unlock its full potential in therapeutic interventions in a variety of pathological process. On this regard, mathematical modelling can help elucidate the mechanisms behind network formation, serving as a proof of concept mean for newly developed theories [99], as well as steer experimental investigations towards the most promising research perspectives.

Mathematical modelling of vasculogenesis
While the mathematical model presented here is the first one formulated to study cluster-based vasculogenesis, many have been proposed to investigate in vitro single-cell vasculogenesis. Most models in the literature refer back to the work of Serini et al. [98] who investigate single-cell vasculogenesis in matrigel. During the early  [16] to study EPC cluster formation in hypoxic gradients (bottom) and in control nonhypoxic environments (top). (B) Corresponding experimental images displaying cluster formation after 48hours in hypoxic environments (bottom), while no clusters formed in nonhypoxic conditions (top). Illustrations and images are taken from Figure 1 in [16] under Creative Commons licence https://creativecommons.org/licenses/by-nc/4.0/. stages of this process, ECs were observed to undergo rapid motion, during which they maintained a round shape, until collision with their closest neighbours (3-6 hours). They then proceeded to reorganise into a continuous multicellular network, which can be represented as a collection of nodes connected by capillary chords (of characteristic length ≈ 200µm), exhibiting a more elongated morphology, until network stabilisation (9-15 hours). Note that cell proliferation and death of mature ECs can be neglected on such a short timescale in vitro. We first review key single-cell vasculogenesis models in the literature, referring the interested reader to more detailed reviews [3,17,31,95], and then discuss model requirements to study cluster-based vasculogenesis.
Partial differential equation models of single-cell vasculogenesis Continuum deterministic models in the extant literature can be categorised as Persistence and Endogenous Chemotaxis (PEC) models [4,29,38,58,62,98], modelling the early stages of single-cell vasculogenesis, and mechanochemical models [64,73,74,76], better suited to describe the late stages of the process, with the exception of the work by Tosin et al. [111] who proposed a comprehensive model.
First introduced by Serini and coworkers [38,98], PEC models consist of a conservation equation for the EC density, a momentum equation describing changes in the EC velocity and a reaction-diffusion equation describing VEGF dynamics. In particular ECs are assumed to undergo persistent motion i.e. movement characterised by inertia in the velocity field, and this movement is modulated by friction between cells and substrate, densitydependent pressure and chemotaxis (with the VEGF being secreted by the cells themselves).
Such models have so far helped to extrapolate information underlying the origin and structure of newly formed vascular networks, for instance relating the characteristic chord length to the VEGF diffusion coefficient D c and decay rate λ c with the formula = D c (λ c ) −1 [4,3,38,98], or investigating the minimum -and maximum -initial cell density required for network assembly [3,29,38,58,98]. In particular, as the initial cell density increases we observe a transition from a phase in which several disconnected structures appear to one in which a single connected structure appears. Detailed analyses of this phenomenon, known as percolative transition, revealed the critical EC density for transition to be around 100 cells/mm 2 [29,38]. Further increasing the initial cell density, instead of a well-defined network with thin chords, one may observe a cell monolayer with void areas called lacunae, a phenomenon referred to as "Swiss-cheese" transition [3,98]. Linear stability analysis revealed chemotaxis to be the key destabilising force for lacunae fomation and pressure to be the main stabilising one [58]. Thanks to the pressure term, the Swiss-cheese transition can be predicted by PEC models about the critical initial cell density of 300 cells/mm 2 [3].
On the other hand, mechanochemical models of vasculogenesis placed special emphasis on the role played by cell traction forces, thus being particularly suited for the later stages of vasculogenesis during which the mechanical interaction between the ECs and the ECM cannot be neglected. Mechanical models or vasculogenesis, following the work of Murray and coworkers [75], consist of a conservation equation for the EC density and one for the ECM density, a force-balance equation to describe changes in the cell-matrix displacement as cells pull on the ECM and, in the case of mechanochemical models, an additional reaction-diffusion equation describing VEGF dynamics [3]. Both ECs and ECM are assumed to be advected according to changes in the cell-matrix displacement, and ECs are also assumed to undergo strain-dependent movement [64,74,76], haptotaxis [76] and chemotaxis [64]. The ECM is modelled as a linear elastic or viscoelastic material, the substratum is introduced via external viscous drag or elastic forces, and cell traction either saturates or goes negative for large densities, the latter understood as modelling density-dependent pressure [76]. Linear stability analysis predicts instability of the homogeneous steady state, and in general lacunae formation, if the cell traction coefficient is high enough or the ECM Young modulus is low enough, and if the initial cell density is below a threshold value [3]. Haptotaxis was also indicated to have a destabilising effect [76], while the restoring forces and long range effects have a stabilising effect [3]. Note that while these models may predict cluster formation in the case of very high cell traction or very low ECM stiffness, such modelling frameworks are not compatible with the application here considered as they do not include any of the key mechanisms of early-stages cluster-based vasculogenesis.
Finally, Tosin et al. [111] proposed a PEC model with an external drag force in the momentum equation, slowing cells down as they move on the ECM, itself modelled as a linear elastic material. The model indicated cell adhesiveness to the substratuum to have a key role in network formation, with too high adhesiveness resulting in no network and too low adhesiveness affecting the stability of chords so that cells would eventually clusterise due to chemotaxis.
Cellular Potts models of single-cell vasculogenesis Despite their lower analytical tractability, cellular Potts (CP) models have been widely used to study single-cell vasculogenesis over the years [17,95]. These are lattice-based models, with simulations obtained via a Monte-Carlo method following energy minimisation principles, which track cell and ECM dynamics at the mesoscale. They are therefore particularly suited to investigate mechanisms occurring at the cell level, such as cell shape and cell-to-cell adhesion, which are difficult to describe with continuum models of macroscale dynamics and which may be crucial for single-cell vasculogenesis as the process involves a low number of cells [67]. These models predicted that elongated, adhesive cells can self-organised into vascular structures [81] and that these structures are faster achieved and stabilised thanks to endogenous chemotaxis [67]. Furthermore, these models confirm a correlation between the characteristic chord length and the VEGF diffusion rate, as previously indicated by PEC models. In the absence of cell elongation, the network formation is compromised and cells form islands. This observation is independent of whether elongation is imposed as a constraint [67] or it is a consequence of strong cell-to-cell adhesion [69], of preferential attraction of cells to elongated structures [106,108,107], or of cell traction, ECM stiffening and durotaxis [86,116]. In the latter case, network formation was predicted for high enough cell traction compared to the ECM Young modulus, just like in mechanochemical models. Alternatively to cell elongation, cell polarisation resulting from cell adhesion-mediated saturation of chemotaxis may prevent clusters from forming [68,70,103]. Overall the work of Merks and coworkers pointed towards cell elongation, or polarisation and endogenous chemotaxis being essential for network formation.
Modelling the early stages cluster-based vasculogenesis While much work has been conducted in order to better understand the mechanisms at the basis of single-cell vasculogenesis, and many of which are relevant to cluster-based vasculogenesis, recent experimental observations [16] seem to suggest that also other processes need to be considered when studying EPC cluster formation during the early stages of cluster-based vasculogenesis. While CP models may be preferred when studying a system comprising of few cells, EPCs actively proliferate and we therefore choose to study this process through the lenses of a continuum deterministic model, which is also more ameanable to analytical investigation. The model should include endogenous chemotaxis, MMPmediated ECM degradation, cell-to-cell adhesion, EPC diffusion, proliferation and death, together with annexed ECM, VEGF and MMP dynamics. Note that while cell proliferation and ECM degradation have been included in previous models of vascular network formation [17,33,50,64,112], these have been formulated to study angiogenesis rather than vasculogenesis, which for instance require different initial conditions and exogenous rather than endogenous chemotaxis.
PEC models included chemotaxis in the momentum equation, but this implies that cells accelerate in chemical gradients, an assumption that might be unrealistic given the highly viscous, non-inertial environment of the ECM [67,70]. We thus choose to follow standard modelling of chemotaxis in the flux term of the EPC density mass-balance equation, as done in mechanochemical models [3,64], following the classic Patlak-Keller-Segel (PKS) model of diffusion and chemotaxis [53,77]. As the PKS model is well-known for having solutions that blow up in finite time, a number of variations have been proposed in the literature [20,57,77], and we will hereby consider a modified version of the forms proposed by Hillen and Painter [47,80] whereby chemotaxis is saturated in a tightly packed environment. While this results in a local description of chemotaxis, various nonlocal ones have also been proposed over the years, reviewed for instance in [26,48,79]. While CP models allowed for a mesoscopic description of cell-to-cell adhesion, during the development of PEC and mechanochemical models this biological process had not yet been explicitly modelled at the macroscopic scale. In 2006 Armstrong and coworkers published a seminal paper [7] proposing a continuum nonlocal description of intrapopulation and interpopulation cell adhesion. This was later adapted by Gerisch and Chaplain [42] to account for cell-to-cell and cell-to-matrix adhesion. The model relies on the presence of an advective flux of the cell density at position x, calculated nonlocally by scouting for availability of adhesion sites in a sensing region centered in x. We refer the interested reader to [46,49,100] for results on the existence and uniqueness of solutions to the model equation proposed by Armstrong et al. [7]. This has become a popular continuum description of cell adhesion, with a number of variations proposed to study tumour invasion and cell movement across the ECM [15,22,25,37,43,78,100], development and cell-sorting dynamics [8,23,43,79]. Attention has also been given over the years to variations of the nonlocal modelling of adhesion by considering different ways to enforce limits on the cell density and furthermore to include nonlinear cross diffusion, that is substituting linear diffusion modelling random movement with nonlinear diffusion modelling movement of cells down a density-dependent pressure gradient, into the overall model, see [21,23,63,72]. The principal effect of these modifications is that they lead to sharper interfaces between the densities of different cell types or the ECM or can even lead to strict segregation of them. Since such effects are not part of the primary dynamics which we will focus on in our model of the early stages of cluster-based vasculogenesis, we will not include these variations here, see however the research perspectives in Section 4.

Study aims and paper structure
The main aim of this work is to develop a continuum deterministic model of vasculogenesis that is suitable for in silico studies of the EPC cluster formation. We thus extend mathematical models of vascular network assembly to investigate the early stages of cluster-based vasculogenesis, during which mechanical interactions between EPCs and the ECM can be neglected. We present a continuum deterministic model of EPC cluster formation, comprising of a system of partial differential equation modelling dynamics such as endogenous chemotaxis, MMP-mediated ECM degradation, nonlocal cell-to-cell and cell-to-matrix adhesion, EPC diffusion, proliferation and death, together with annexed ECM, VEGF and MMP dynamics. Our model provides a theoretical basis for a comprehensive description of the mechanisms underlying cluster formation. In such perspective, another objective of this study is to clarify the role played by different dynamics and elucidate the determinants of cluster size with the introduction of appropriate metrics. We investigate this in primis by means of a linear stability analysis and numerical simulations relying on a baseline parameter set drawn from the literature. We then proceed to provide a first overview of potential model behaviour and an indication of the importance of various parameters or processes on the pattern formation potential. We do this by simulating the model using the baseline parameter set but with different parameters, one at a time, taking values from a suitable range of values as mostly identified from existing literature. Finally, in view of our findings, we seek to identify crucial modelling and experimental research perspectives, together with mathematical investigations that could yield significant results and complement experimental studies of cluster-based vasculogenesis.
In Section 2 we describe the mathematical model in detail, together with the analytical method used to investigate the system and the numerical method used to obtain solutions, and introduce the metrics defined to study cluster size. In Section 3 we report the analytical results, followed by detailed numerical investigations of the cluster formation process and the size of clusters for the one-dimensional (1D) and two-dimensional (2D) problems, which are qualitatively compared with the experimental findings of Blatchley et al. [16]. In Section 4 we summarise the main findings and discuss future research perspectives.

Mathematical model
Let t ∈ R ≥0 indicate time and x ∈ R 2 position in space. We will consider the dynamics in a 2D spatial domain Ω ⊂ R 2 , as well as the corresponding 1D problem, for which we make use of the notation x ∈ R to indicate space. Unless indicated otherwise, all definitions introduced in this section for the 2D problem hold for the corresponding 1D one. The density of EPCs at time t and position x is given by n(t, x), in units of cell cm −3 , and the ECM density by ρ(t, x), in units of nM. Similarly we indicate the concentration of a matrix-degrading enzyme, such as MMP-1, by m(t, x), in units of µg cm −3 , and that of a chemotactic agent, such as VEGF-A, by c(t, x), in units of ng cm −3 . We also introduce the vector of dependent variables v(t, x) := n(t, x), ρ(t, x), m(t, x), c(t, x) . Our model then consists of a system of mass-balance equations, one for each dependent variable introduced.

EPC density dynamics
The mass-balance equation for the density of EPCs is of the form: where J d (n) models spatial diffusion to account for random movement of cells, J c (n, c) indicates the chemotactic flux in response to VEGF gradients and J a (v) the advective flux due to cell-to-cell and cell-to-matrix adhesion, while the term on the right-hand side of the equation models cell proliferation and death. In particular, we consider a modified version of the standard logistic growth, as proposed by Gerisch and Chaplain [42], in which EPCs proliferate at rate p ∈ R ≥0 and die due to competition for space -occupied by both cells and ECMand resources. Parameters ϑ 1 ∈ R >0 and ϑ 2 ∈ R >0 indicate the fraction of one unit volume of physical space occupied by EPCs at unit density and by ECM at unit density respectively, such that (ϑ 1 n + ϑ 2 ρ) indicates the total fraction of locally occupied space.
Spatial diffusion and chemotaxis Following the motivations introduced in Section 1.2, we make use of the following definitions for the diffusive and chemotactic flux terms: In definitions (2) spatial diffusion follows Fick's law, with diffusivity D n ∈ R ≥0 , while chemotaxis is modelled as an advective flux of cells up the gradient of the chemotactic agent concentration c, modulated by the chemotactic sensitivity of cells. This is proportional to the chemotactic sensitivity coefficient χ ∈ R ≥0 and the function f (n, ρ), which accounts for volume exclusion. This is given by where we have used (·) + := max(0, ·) -see also [42,47,80] . According to definition (3), the chemotactic sensitivity of the cells is proportional to the locally available space. If the local space is overcrowded, cells will struggle to sense the chemotactic gradient, to the extent that if the space is locally full no chemotaxis can occur.
Cell-to-cell and cell-to-matrix adhesion Following the motivations introduced in Section 1.2, we model cell-to-cell and cell-to-matrix adhesion in the continuous form proposed by Gerisch and Chaplain [42]. The advective flux of the EPC density n due to cell-to-cell and cell-to-matrix adhesion in (1) is therefore given by where A[v(t, ·)] is the adhesion velocity at some point x, that is the velocity of cells at x due to adhesive interactions with their environment, which is an operator acting on v(t, ·) defined as a function of x. For problems in 1D this is defined as while for problems in 2D it is defined as In equation (5) η(j) = (−1) j with j = 0, 1 indicates the 1D right and left unit outer normal vector, while in equation (6) the 2D unit outer normal vector corresponding to angle θ is given by η(θ) = (cos θ, sin θ) . In definitions (5) and (6), we have that the sensing region of cells at position x ∈ R d (d = 1, 2) is the d-dimensional ball centred in x with radius R > 0, called the sensing radius. Then Γ(r) is the radial dependency function, indicating how strong the adhesion velocity at a point x is influenced by points at a distance r ≤ R from the centre x of the sensing region. Since this should not alter the magnitude itself of the adhesion velocity, Γ(r) must satisfy in the 1D and 2D problems respectively, i.e. in equations (5) and (6) respectively. Assuming Γ(r) decays linearly with r to be zero at the boundary of the sensing region, we here make use of the forms proposed by Gerisch and Chaplain [42] in the 1D and 2D problems respectively, chosen so that (7) is satisfied. Finally, in both (5) and (6) the term g(v(t, ·)) represents the nonlocal impact of the system's state at some point within the sensing region at x on the velocity of the cells at x due to adhesions to other cells or the ECM in the sensing region at x. This is given by in which S nn ∈ R ≥0 and S nρ ∈ R ≥0 are the cell-to-cell and cell-to-matrix adhesion coefficients respectively, while the function f (n, ρ) is defined in (3). Under definition (9), the velocity of cells at position x in the direction of a point -say -y in the sensing region of x due to cell-to-cell adhesion is directly proportional to the cell density n at y, and that due to cell-to-matrix adhesion is directly proportional to the ECM density ρ at y. This is because a higher cell or ECM density correlates with a higher number of adhesion sites. Meanwhile, under definitions (9) and (3), the velocity is also proportional to the available space at position y. This accounts for volume exclusion, as cells will be unable to sense adhesive sites at spatial locations with high densities of EPCs and/or ECM and, hence, will not migrate in those directions.

ECM density dynamics
We let the ECM be degraded by matrix-degrading enzymes at a rate γ ∈ R ≥0 and account for ECM remodelling at a rate µ ∈ R ≥0 , resulting in the following mass-balance equation for the ECM density: where the parameters ϑ 1 and ϑ 2 have already been introduced in Section 2.1.1, indicating ECM remodelling is here understood as a restructuring phenomenon that only occurs if space is available [37]. This is assumed to be independent of the cell density n as ECM remodelling in vivo is generally mediated by other cell types present in the tissue, such as mesenchymal stem cells, fibrocytes and fibroblast [12,14,36,66,104], which we do not include in our modelling framework.

MMP concentration dynamics
We let the matrix-degrading enzyme (MMP) be produced by the EPCs at a rate α m ∈ R ≥0 , undergo Fickian diffusion with diffusivity D m ∈ R ≥0 , and decay at rate λ m ∈ R ≥0 . Then the MMP concentration m(t, x) satisfies:

VEGF concentration dynamics
Similarly to the matrix-degrading enzyme, we let the chemotactic agent (VEGF) be produced by the EPCs at a rate α c ∈ R ≥0 , undergo Fickian diffusion with diffusivity D c ∈ R ≥0 , and decay at rate λ c ∈ R ≥0 . This results in the following mass-balance equation for the VEGF concentration c(t, x):

Boundary conditions
While equation (10) describes the dynamics of the ECM in the closed spatial domainΩ = Ω ∪ ∂Ω, equations (1), (11) and (12) are posed on the open set Ω, and are complemented with zero-flux boundary conditions. These boundary conditions imply that no mass is exchanged with the outside of the spatial domain, i.e. we have a closed system. For the nonlocal terms (5) and (6) this means that they cannot sense the system's state outside of the spatial domain Ω. We thus impose that in these terms the function g(v(t, ·)) equals zero if it is to be evaluated for a point outside of Ω, see also [37], and thus obtain well-defined nonlocal terms (5) and (6) throughout the spatial domain. The above definition of the nonlocal terms near the spatial domain boundary resembles what is termed repulsive (or naive) boundary condition in [46]. This reference also gives further possibilities for the definition of nonlocal terms on bounded domains and discusses their mathematical and biological implications.

Initial conditions
As proposed by Serini et al. [98], we construct the initial conditions to mimic sparsely distributed cells on the ECM. In particular, the initial cell density is given by the sum of K ∈ N randomly distributed bell-shaped bumps. In particular, we construct these bumps as Gaussian-like functions with maximum height and full width at half maximum (FWHM) both equal to the size of an average cell diameter a ∈ R >0 . We let the initial ECM density be constant and the initial MMP and VEGF concentrations be null. Thus we have where G i (x; a) indicates the Gaussian-like function centered at the (randomly selected) x (i) ∈ Ω and is given by Initial conditions (13) and (14) are defined for the 2D problem and we make use of analogous definitions for the 1D problem, obtained by simply replacing x ∈ Ω ⊂ R 2 with x ∈ Ω ⊂ R. In equation (14) we have used the formula FWHM= 2 √ 2 ln 2σ where σ is the standard deviation of the Gaussian.

Nondimensional model and baseline parameter set
We nondimensionalise the system of equations (1) and (10)- (12), together with definitions (2)-(9), (13) and (14), by lettingt We use L = 0.1 cm as characteristic length scale, in accordance with previous vasculogenesis works [64,98] and for easy visual comparison with the experimental results reported by Blatchley and coworkers [16]. We then take reference time scale τ : and we take ϑ 1 = 10 −9 cm 3 /cell, the average volume occupied by an endothelial cell [89]. We use a reference ECM density of P = 10 −1 nM [5,6,110] and define the parameter ϑ 2 := P −1 . We take the reference VEGF density to be C = 20 ng cm −3 , in the range of values generally considered in in vitro set ups [45,60,98]. Finally, Blatchley et al. [16] reported concentrations of MMP-1 in the range 1 − 100 µg ml −1 , so we take the intermediate concentration as reference MMP density, i.e. M = 10 µg cm −3 . Let us introduce the following nondimensional parameters: Then the overall nondimensionalised system becomes, dropping hats for convenience, where ( (6) in 2D, with Γ(r) still defined as in (8), and g(v) is given as in (9) with f (n, ρ) now given by The parameters R, S nn and S nρ are the nondimensional ones introduced above, and the system (15) is complemented with initial conditions (13)- (14), in which the parameters a and ρ 0 now corresponds to the nondimensional ones introduced above. The baseline parameter set, with the corresponding nondimensional (ND) parameter values, is reported in Table 1 -see A for details.

Linear stability analysis
We perform a linear stability analysis (LSA) on the spatially homogeneous steady states, say,v of the nondimensional model (15)-(16), to gain insights into possible mechanisms responsible for aggregation dynamics. During the LSA we first introduce a small spatially homogeneous perturbation v =v +ṽ(t), with |ṽ| 1, in (15) and linearise. By assuming the small perturbation is proportional to exp (σt), we derive a dispersion relation for σ and impose Re(σ) < 0, to ensure spatially homogeneous steady states are stable to spatially homogeneous perturbations. We then repeat these steps under a small spatially inhomogeneous perturbation v =v +ṽ(t, x), with |ṽ| 1, assuming it is proportional to exp (σt + ikx) -orṽ(t, x) ∝ exp (σt + i k · x) in 2D. Once we obtain a dispersion relation σ(k 2 ) -where k 2 = |k| 2 in 2D -we study the conditions under which Re(σ(k 2 )) > 0 for some k 2 , as in such regimes we expect spatially inhomogeneous perturbations to grow in time and patterns to arise. We conduct this analysis both in 1D and 2D for the model (15)- (16), as well as for the corresponding problem in absence of volume exclusion, i.e. substituting definition (16)

Numerical method
We solve the nondimensional system (15)- (16) in Ω = (0, 1) in 1D and Ω = (0, 1) × (0, 1) in 2D, with zeroflux boundary conditions and initial conditions (13)- (14). All numerical simulations have been performed in Matlab. The numerical scheme follows the method of lines by first discretising the nonlocal model in space (with 1000 grid cells in 1D and 100 × 100 grid cells in 2D), yielding an initial value problem for a large system of stiff ordinary differential equations. This system is then solved using the linearly-implicit time integration scheme ROWMAP [117], implemented in a Fortran subroutine and called from Matlab. For the discretisation in space we use a second-order finite volume approach which makes use of flux-limiting in order to ensure an accurate and at the same time non-negative approximation of the taxis and adhesion terms; note that here in particular we employ the Koren flux limiter (k = 1/3, δ = 0.25). This finite-volume discretization in space, except for the approximation of the nonlocal term, is fully described in [41] and references therein. A detailed account of the approximation of the nonlocal term in a periodic boundary condition setting is given in [40]. The essential observation there is that, for the simultaneous evaluation of the nonlocal term on all grid cell interfaces, a matrix-vector product with a circulant matrix (spatially 1D case) or a block-circulant matrix (spatially 2D case) must be evaluated and that this can be done efficiently using Fast Fourier Transform (FFT) techniques. In fact, this efficient approximation of the nonlocal term, simultaneously on the full computational grid, is key to an overall efficient numerical scheme for the full PDE system. In our case here, however, we have that the nonlocal term comes with zero-flux boundary conditions and the matrix-vector product in this case involves matrices with Toeplitz structure -see also [37]. In this case, FFT techniques cannot be applied directly. However, these (block) Toeplitz matrices have a banded structure with a band width which is small thanks to the fact that the size of the sensing region of the nonlocal term is small compared to the size of the domain Ω. Thus these (block) Toeplitz matrices can be embedded in slightly larger (block) circulant matrices and the desired result of the matrix-vector product can again be obtained efficiently using FFT techniques. We elaborate on this embedding in the Supplementary Material. Recent work by Colombi et al. [28] considers a nonlocal term similar to ours and investigates numerical approaches to its approximation. Note, however, that in their case the nonlocal term is evaluated in a single spatial location (changing with time) as opposed to the computationally much more involved task of approximating the nonlocal term on all grid cell interfaces simultaneously as in our case.

Cluster size metrics
In order to gain insight into the role played by different biological, chemical and mechanical factors in dictating the cluster size, we define two different measures of cluster size with complementary information, as similarly done in [82]. We define these for the nondimensional 1D problem, but analogous definitions can be considered for the 2D problem. Assume that Q ∈ N clusters have formed at time t = T and let ω ⊂ Ω be the subdomain supporting these clusters, i.e.
Then ω can be partitioned into Q subdomains ω 1 , ..., ω Q , i.e. we have where each ω i (i = 1, ..., Q) corresponds to the support of a cluster. We let the average cluster width W and average cluster compactness C be defined by Under definitions (18)-(20), the width W i (i = 1, ..., Q) of each cluster is a measure of the length of its support, which may be understood as an indicator of the diameter of the cluster assuming the 1D case is reflective of the corresponding 2D problem. Note that the analogous 2D definition to (20) would inform us on the area covered by each cluster, from which the average cluster diameter could be calculated. However, this would need to be complemented with an additional metric for cluster elongation (e.g. the ratio between the diameter of the circle circumscribing ω i and that of the one inscribed in ω i ), in order to obtain an exhaustive description of the cluster's topology. Under definition (21), in which n(T, x) is the nondimensional cell density, the compactness C i (i = 1, ..., Q) of each cluster is a measure of the average cell density within cluster i. Cluster compactness allows us to distinguish between simple cell aggregates and well-defined clusters, identified as such only if C is higher than 0.5, corresponding to at least half the local volume being occupied by cells. Under the choices of nondimensionalisation and of spatial domain for the numerical simulations, we expect 0 ≤ W i ≤ 1 andunder cell incompressibility assumptions -0 ≤ C i ≤ 1 for all i = 1, .., Q. Figure 3 summarises the biological interpretation of possible combinations of W and C. In practice, the cluster domains ω i used to calculate W and C are identified via image segmentation by thresholding, with a threshold value for the cell density set to be 10 −4 (two orders of magnitude smaller than the nondimensional average cell diameter a), below which n(t, x) is approximated to zero.

Results
As outlined in Section 1.3, we here investigate the determinants of cluster formation and cluster size, by means of a linear stability analysis and numerical simulations, thank to which we provide an overview potential model behaviours under one-at-a-time perturbations from the baseline parameter set. We report in Section 3.1 the key results of the LSA, as described in Section 2.3, and in Sections 3.2-3.5 the results of the numerical simulations obtained using the method described in Section 2.4. In Section 3.2 we report the results of the 1D model under the baseline parameter set in Table 1 and qualitatively investigate the determinants of cluster formation in Section 3.3. In Section 3.4 we report the results of the parametric analysis conducted to elucidate the relation between different parameters and the average cluster width and compactness, as defined in Section 2.5. Finally, we report the results of the 2D model in Section 3.5.

Linear stability analysis results
We first conducted a LSA on the biologically significant spatially homogeneous steady states of the system (15) under definitions (5)- (9) and (16), both in 1D and 2D. Then, we investigated any changes in the LSA results in the absence of volume exclusion in the chemotactic and cell adhesion terms of equation (15) 1 , that is under definition (17) for f (n, ρ) instead of (16). The full LSA can be found in Section S1 of the supplementary document 'SuppInfo' with the key results being as follows. We find that the two spatially homogeneous steady states are such that either the whole domain is solely occupied by cells (n = 1 andρ = 0) or solely occupied by ECM (n = 0 andρ = 1), the latter of which we do not consider further as it is not biologically relevant for our current study. In absence of ECM degradation (γ = 0), in addition to the two above, there are infinitely many "intermediate" spatially homogeneous steady states satisfying 0 <n,ρ < 1, withn +ρ = 1. These "intermediate" steady states are linearly stable to small spatially inhomogeneous perturbations. When we neglect the influence of volume exclusion in the chemotactic and cell adhesion terms of equation (15) 1 , we find that chemotaxis and cell-to-cell adhesion can combine to make the biologically significant steady states (n > 0), which are the same as in the presence of volume exclusion, unstable to small spatially inhomogeneous perturbations if the cell density is sufficiently large (see the 'SuppInfo' document for more details). Thus we expect patterns to arise when cell-to-cell adhesions and chemotactic mechanisms dominate the dynamics, which may occur, for instance, when the domain is not densely packed with cells and ECM. We therefore take the initial ECM density in (13) to be ρ 0 = 0.5 in our nondimensional baseline parameter set, in order to increase the chances that clusters form in the model.

Cluster formation under the baseline parameter set
We report in Figure 4 the cell density n(t, x) obtained from numerical simulations of the 1D model under the baseline parameter set and in Figure 5 the corresponding cluster width W and compactness C. The plots displayed in these figures indicate that our model predicts three stages of cluster formation in 1D. First EPCs form aggregates which reach minimum cluster width of about W = 0.2 around t = 17 (cf. left panel in Figure 4 and Figure 5). Then the cells in these aggregates continue to proliferate increasing their compactness while keeping the cluster width unchanged up to about t = 50 (cf. central panel in Figure 4 and Figure 5). Finally, the cells continue to proliferate until the whole domain is occupied by cells (cf. right panel in Figure 4 and Figure 5). The simulated dynamics of cluster formation nicely match the experimentally observed ones by Blatchley et al. [16], see Figure 1, although in our simulations they occur on a slightly slower timescale. In fact in [16] the Figure 4: Plots of the cell density n(t, x) obtained solving the system (15), together with definitions (5), (8) and (16), initial conditions (13) and (14)  minimum cluster size was reached around 24 hours (about nondimensional time t = 8.7) and kept unchanged while clusters increased compactness up to 48 hours (t = 17.3) before late stage dynamics kick in. The numerically obtained cluster width corresponds to about 200 µm, which is within the range of cluster diameters observed in [16]. Finally, after t = 50 we observe the equivalent of whole tissue invasion, suggesting our modelling framework cannot properly capture cluster stabilisation in the long run, at least in the 1D case.  Table 1, reported in

Qualitative investigation of the determinants of cluster formation
In view of the results reported in Section 3.2, we now investigate the role played by chemotaxis, ECM degradation and cell-to-cell adhesion in cluster formation by varying the relevant parameter values and observing changes in the solution up until t = 50, starting from the same initial conditions considered in the previous section. The long-time dynamics (up to t = 400) of the solutions plot in Figure 6 are reported in the supplementary Figure S2.
The primary role of chemotaxis and ECM degradation The plots reported in the second row of Figure 6 reveal the role of chemotaxis and ECM degradation in cluster formation, according to the dynamics described in Section 2.1. Under the baseline parameter set we observe cluster formation (cf. first plot in second row of Figure 6), as discussed in the previous paragraph and summarised in Figure 4. In absence of ECM degradation, even though cell aggregates of the same width form, the maximum cell density remains below 0.3 (cf. second plot in second row of Figure 6) and actually decays over longer periods of time (see supplementary Figure S2). On the other hand, in absence of chemotaxis no cell aggregation occurs and we either observe total invasion of the domain by the cells (cf. third plot in second row of Figure 6) or, in absence of ECM degradation, a simple spatial redistribution of the cells (cf. fourth plot in second row of Figure 6). These results indicate that chemotaxis and ECM degradation are both crucial to cluster formation, with chemotaxis playing a key role in cell aggregation and ECM degradation being responsible for these aggregates to grow into well-defined and compact clusters.
The secondary role of cell-to-cell adhesion Let us now compare the plots in the second row of Figure 6 with those in the rest of the figure, which have been obtained by varying the cell-to-cell adhesion coefficient S nn , the value of which was chosen a priori due to lack of proper estimates in the current literature. We immediately observe that for small values S nn ≤ 0.1, under the baseline parameter set, cell-to-cell adhesion does not play any particular role in cluster formation, as demonstrated by the fact that the simulations in absence of cell-to-cell adhesion closely resemble those with S nn = 0.1 (cf. first and second row of Figure 6 and Figure S3 in 'SuppInfo' document). On the other hand, increasing the order of magnitude of the cell-to-cell adhesion coefficient results in the initial formation of small-scale aggregates. Moreover, the maximum density reached by these aggregates increases as S nn increases (cf. third and fourth row of Figure 6). This however does not seem to affect the long-time dynamics of the solution, which remain analogous to those described in the previous paragraph. For instance, under the baseline parameter set except for S nn = 1, the initial aggregates merge -likely due to chemotaxis -into wider clusters and then proceed to invade the whole tissue (cf. first plot in third row of Figure 6 and of Figure S2). For S nn = 10 we do not yet observe small aggregates merging into larger clusters at t = 50 -likely due to cell-to-cell adhesion overpowering chemotaxis -and tissue invasion is simply delayed (cf. first plot in fourth row of Figure 6 and of Figure S2). Overall these results suggest that the continuum nonlocal description of cell-to-cell adhesion considered in our model, while it may capture the aggregating effect of cell-to-cell adhesion for S nn high enough, does not capture the stabilising effect that we are seeking for in this modelling framework.

Quantitative investigation of the determinants of cluster size
We study changes in the measures W and C, defined according to (20) and (21), at t = 50 under variations of each parameter in the baseline parameter set. In particular, we consider the effect of halving and doubling the magnitude of each parameter in equations (15) 1 and (15) 2 in Figure 7, and those in equations (15) 3 and (15) 4 in Figure 8. Each boxplot had been obtained using data from 100 simulations under the same parameter set, staring from randomised initial conditions as in (13) (15), together with definitions (5), (8) and (16), initial conditions (13) and (14) 1 , complemented with zero-flux boundary conditions, in absence of cell-to-cell adhesion, i.e for S nn = 0: under the baseline parameter set (first column), in absence of ECM degradation, i.e. for γ = 0 (second column), in absence of chemotaxis, i.e. χ = 0 (third column), and in absence of both chemotaxis and ECM degradation, i.e. χ = γ = 0 (fourth column). Second, third and fourth rows: Same as first row but in the presence of cell-to-cell adhesion, with S nn = 0.1 (second row), S nn = 1 (third row) and S nn = 10 (fourth row) respectively. plot). The nondimensional width corresponds to a diameter in the range 150 − 410µm, which is consistent with the experimentally observed cluster sizes in [16] (cf. Figure 1).
The role of chemotaxis Boxplots of the width W and compactness C of clusters for different values of the chemotactic sensitivity χ are displayed in Figure 7b. While higher values of χ correlate with slightly smaller clusters, yet maintaining a mean width around 0.2, lowering the magnitude of χ seems to result in a wider range of values of W with much higher median, as well as higher compactness C. This supports the notion that lowering the chemotactic sensitivity of the cells hinders cluster formation and -assuming all other dynamics are present -fosters tissue invasion, which is in line with the results presented in Section 3.3.
The role of ECM degradation Boxplots of the width W and compactness C of clusters for different values of the ECM degradation rate γ are displayed in Figure 7f. Notice that, while the median width W is maintained around 0.2, increasing the magnitude of γ results in a smaller range of values measured for both W and C, as well as higher values of cluster compactness. In addition, this trend suggests that further decreasing γ will lead to higher values of W and lower values of C measured (see supplementary Figure S4c in 'SuppInfo' document), which is in line with the observed dynamics in absence of ECM degradation in Section 3.3 (cf. Figure 6, second and fourth columns). Overall, this highlights the key role ECM degradation has in cluster formation, establishing a relation between the rate γ and the cluster compactness C.
The role of matrix remodelling We see in Figure 7g that increasing the magnitude of the matrix remodelling rate µ results in a larger range of values measured for both W and C, as well as lower values of cluster compactness, while the median width W is maintained around 0.2. These are opposite effects to those obtained increasing its degradation rate γ (vid. Figure 7f), which is coherent with the opposite nature of these dynamics.
The role of cell-to-cell adhesion The plots in Figure 7c confirm that small variations in the cell-to-cell adhesion coefficient S nn do not influence the size of clusters, as expected from the results in Section 3.3. Additional numerical tests considering different orders of magnitude of S nn (see supplementary Figure S1b in 'SuppInfo' document) revealed a slight increase in the median W , probably due to the initial presence of smaller clusters (cf. Figure 6, third row) which at t = 50 may still be in the process of merging depending on their spatial distribution. These results do not suggest that cell-to-cell adhesion is a key mechanism in determining cluster -and therefore network -topology.
The role of cell diffusion and cell-matrix interactions The boxplots of W and C of clusters for different values of the cell diffusion coefficient D n , displayed in Figure 7a, show that changing the magnitude of D n has very little effect on the cluster compactness C and median width W . There is however an increasing range of values of W measured as D n increases, indicating that low diffusivity correlates with more precisely defined clusters, while high diffusivity results in more variability in cluster size. Moreover, this variability allows for larger values of W to be measured, suggesting that much higher diffusivity may result in tissue invasion (see supplementary Figure S4a in 'SuppInfo' document). The same observations can be conducted on the boxplots in Figure 7d, obtained by varying the magnitude of the cell-to-matrix adhesion coefficient. This is in line with the notion that lower cell-matrix interactions facilitate cluster formation, while much higher cell-matrix interactions promote tissue invasion.
The role of cell proliferation In Figure 7e we see that slower proliferation -i.e. lower p -correlates with a wide range of lower values of cluster compactness C measured, while the median width is maintained around 0.2. On the other hand, faster proliferation -i.e. higher p -results in a wide range of higher values of cluster width W and a small range of high values of compactness C. Note that these data portray different stages of the cluster formation process, as demonstrated in Figure 5 (central panel): initially aggregates form without being very condensed (low-to-medium W and C), then they increase their compactness while keeping steady width (low W and high C), and eventually grow further invading the surrounding space (medium-to-high W and high C). This suggests that the rate of proliferation of EPCs might play a key role in determining the speed of the cluster formation process.
The role of initial ECM density The plots in Figure 7h indicate that changes in the initial ECM density ρ 0 does not affect the long-time spatio-temporal dynamics. In fact, while different values of ρ 0 results in slightly different values of the median width and compactness of clusters observed at t = 50, these values still capture the same biological scenario (i.e. cluster are recognisable at t = 50).
The role of VEGFs dynamics Figure 8e displays changes in cluster width W and compactness C as the rate production of VEGFs α c varies. Note that these results well mirror those obtained by varying chemotactic sensitivity (cf. Figures 7b and 8e), which is coherent with the notion that higher VEGF production rates correlate with stronger chemotactic dynamics, already established to play a key role in cluster formation. In addition, the size of clusters seems to be proportional to the VEGF diffusion coefficient D c , as demonstrated by the plots in Figure 8d in which we see that increasing the magnitude of D c results in higher W and C. This trend, however, suggests that much higher values of D c may result in tissue invasion, rather than cluster formation (see supplementary Figure S4e in 'SuppInfo' document). Finally, as suggested by Figure 8f, changes in the VEGF decay rate λ c do not seem to particularly effect cluster size over the range of parameter values considered. Figure 8b displays changes in cluster width W and compactness C as the rate production of MMPs α m varies. Similarly to what was observed for α c in relation to χ, we see that these boxplots closely resemble those obtained varying γ (cf. Figures 7f and 8b), which is coherent with the notion that higher MMP production rates correlate with stronger degrading dynamics, already established to be responsible for turning aggregates into clusters -that is, increasing their compactness. In addition, we see in Figure 8c that increasing the magnitude of the MMPs decay rate λ m yields opposite effects to those obtained increasing their production α m , further confirming the role MMP-mediated ECM degradation has in cluster formation. On the other hand, changes in the MMPs diffusivity D m do not seem to affect cluster size, as demonstrated in Figure 8a -verified under different orders of magnitude of D m (see supplementary Figure S4d in 'SuppInfo' document). This suggests not much will be gained by distinguishing between membrane-bound and diffusive MMPs, at least in this modelling framework and within the range of parameter values considered.

Two-dimensional clusters
Let us now consider the 2D problem. In this section we focus on the most interesting results obtained in Section 3.4 and investigate the role played by chemotaxis, ECM degradation and cell proliferation in the formation of 2D clusters. We also consider the results in absence of ECM remodelling, discussing its biological interpretation in relation to in vitro and in vivo assays.
2D clusters under the baseline parameter set Under the baseline parameter set reported in Table 1, 2D cluster formation follows slightly different spatio-temporal dynamics to those observed in the 1D case, as demonstrated by the plots in the second row of Figures 9, 10 and 11. At t = 20 the cell density has already reached maximum local compactness in some regions, while aggregation dynamics are still at their early stages (cf. first panel in second row of Figures 9, 10 or 11), and the minimum cluster size is observed no earlier than t = 140 (cf. fourth panel in second row of Figures 9, 10 or 11). The cluster observed at this stage has a nondimensional diameter of about 0.2, which agrees with the results in 1D and is therefore also consistent with experimental observations reported by Blatchley et al. [16]. After the cluster has formed, no tissue invasion is observed -see supplementary Figure S5 in 'SuppInfo' document.
The role of ECM degradation in 2D The plots reported in Figure 9 demonstrate that ECM degradation promotes the formation of 2D clusters, as predicted by the results of the 1D model presented in Sections 3.3 and 3.4. Lower ECM degradation rates correlate with slower aggregation dynamics and lower compactness of      (20) and (21), under variations from the baseline parameter set in Table 1  such aggregates (cf. first row of Figure 9) -to the extent that in absence of ECM degradation no clusters form (see supplementary Figure S6 in 'SuppInfo' document) -and higher ECM degradation rates correlate with faster cluster formation with well-defined and compact clusters observed at much earlier times (cf. third row of Figure 9). As for the 1D model, we obtain the same results under analogous variations of the MMP secretion rate α m (see supplementary Figure S7 in 'SuppInfo' document). The role of chemotaxis in 2D The plots reported in Figure 10 demonstrate that endogenous chemotaxis promotes aggregation dynamics, as predicted by the results of the 1D model presented in Sections 3.3 and 3.4, and reveal that chemotaxis is an important determinant of 2D cluster topology. Lower values of the chemotactic sensitivity χ correlate with slower and weaker aggregation dynamics, so that no well-defined clusters can be observed (cf. first row of Figure 10). Higher values of χ correlate with faster and stronger aggregation dynamics, with well-defined clusters observed at earlier times (cf. third row of Figure 10) and cluster diameter -once the clusters have reached minimum size -smaller than that observed with lower values of χ (cf. last panel in the second row and last panel in the third row of Figure 10). As for the 1D model, we obtain the same results under analogous variations of the VEGF secretion rate α c (see supplementary Figure S8 in 'SuppInfo' document).
The role of cell proliferation in 2D While the results of the 1D model presented in Section 3.4 seemed to suggest that the speed of the cluster formation process is proportional to the cell proliferation rate p, the plots reported in Figure 11 in a 2D framework indicate otherwise. While the rate of cell proliferation may influence how the cells respond to spatial gradients at the beginning of the cluster formation process (cf. plots in the second column of Figure 11), it does not seem to affect the overall spatio-temporal dynamics of 2D cluster formation. Upon these considerations, the results displayed in Figure 7e may simply be a 1D projection of the 2D dynamics occurring around the same time -consider for instance a 1D cross section (e.g. x 1 = 0.5) of the plots in the second column of Figure 11.
The role of matrix remodelling and its relation to in vitro studies While matrix remodelling naturally occurs in vivo thanks to the presence of other cells in biological tissue, this is not generally observed in in vitro assays, which occur in isolated environments. Thus, in order to compare our results with those of in vitro studies, we chose to set the ECM remodelling rate µ = 0, the results are presented in Figure 12. We observe that clusters still form in the absence of ECM remodelling, however they are smaller than in the baseline parameter set (cf. last panel in the second row of Figures 9-11 and that of Figure 12). The smaller cluster sizes may result from the more dominant role played by cell-to-cell adhesion following ECM degradation: in absence of ECM remodelling, volume exclusion is less likely to affect the cell dynamics and, at the cell boundaries, cell-matrix interactions become negligible compared to cell-cell interactions. We would then expect to observe clusters of a smaller diameter but higher maximum density, but because of death due to competition for space we instead observe a loss of cell mass. While these results are mathematically consistent with our model set up, the observed cluster size under the baseline parameter set was more coherent with in vitro experimental observations in [16] than that observed for µ = 0. This could be explained by considering that our 2D set up is meant to reflect dynamics occurring in a 2D horizontal cross section of the 3D experimental domain reported in [16], illustrated in Figure 1A. As the cells degraded the matrix and reorganised into clusters, they also fell towards the bottom of the hydrogel, as illustrated by the experimental data reported in Figure 13, in regions where it was not yet degraded. Thus, following the horizontal plane intersecting the cluster's center, we would likely observe section crossing the center of a cluster, we would indeed observe cell-independent matrix remodelling.

Discussion and research perspectives
Despite the great progress made in the past 20 years in understanding the mechanisms behind EPC cluster-based vasculogenesis, much more needs to be achieved in order to unlock its full therapeutic potential. Mathematical modelling provides theoretical means to shed light on the otherwise hidden role played by underlying dynamics in the origin and structure of the emergent vascular network, as previously achieved in the study of mature EC vascular network formation (single-cell vasculogenesis). We therefore formulated a nonlocal PDE model of EPC cluster formation during the early stages of vasculogenesis, including mechanisms such as ECM degradation, cell proliferation and cell-to-cell adhesion, which were recently found to distinguish cluster-based vasculogenesis Figure 12: Plots of the cell density n(t, x) obtained solving the system (15), together with definitions (6), (8) and (16), initial conditions (13) and (14)   from single-cell vasculogenesis [16]. Thanks to the introduction of appropriate metrics of cluster width and compactness, we investigated the role played by underlying dynamics in facilitating cluster formation, regulating the speed of the cluster-formation process and the size of clusters-see Table 2 for a summary of the results of the parametric analysis of the 1D problem. Furthermore, we verified that the key observations from the parametric analysis for the 1D model still hold in a 2D framework, with the exception of the role played by EPC proliferation. Parameter Main effect of increasing the parameter value on cluster size Ref. Figure  D n Little effect favouring tissue invasion 7a, S4a χ Better defined clusters (smaller range of width W measured) 7b S nn Little to no effect 7c, S4b S nρ Little effect favouring tissue invasion 7d p Increased cluster formation speed (refuted in 2D) 7e γ Better defined clusters (higher median compactness C measured) 7f, S4c µ Less well defined clusters (lower median compactness C measured) 7g ρ 0 Little to no effect 7h D m Little to no effect 8a, S4d α m Better defined clusters (higher median compactness C measured) 8b λ m Less well defined clusters (lower median compactness C measured) 8c D c Bigger clusters (higher median width W and compactness C measured) 8d, S4e α c Better defined clusters (smaller range of width W measured) 8e λ c Little to no effect 8f

Chemotaxis, degradation and their link to hypoxia
Our results confirmed the role played by ECM degradation in the formation of EPC clusters both with the 1D and 2D models, providing additional theoretical support to the mechanism of EPC cluster formation proposed by Blatchley et al. [16]. For example, the in vitro experiments indicated that no clusters form in the absence of MMPs, which is precisely what we discover from the numerical solutions of our model. In addition, the investigation conducted in Sections 3.4 and 3.5 indicates that the speed of cluster formation is proportional to the rate of ECM degradation -or the rate of MMP secretion by the cells -which nicely agrees their experimental observations. Our numerical results further highlighted that ECM degradation alone may not suffice to explain the formation of clusters, as endogenous chemotaxis was shown to be responsible for aggregation dynamics, without which the cells would simply invade the whole tissue. In addition, the investigation conducted in Sections 3.4 and 3.5 suggests that the size of clusters is (inversely) related to the chemotactic sensitivity of the cells -and the rate at which they secrete VEGF -indicating chemotaxis may be a key determinant of cluster topology. Note that, in view of the experimental evidence presented in [16] and [1], the MMP production rate α m and the VEGF production rate α c may be correlated with the local level of hypoxia. It would therefore be interesting to let these production rates be given as functions of the local oxygen concentration. The results here obtained suggest that our model would then predict cluster formation to be fostered by higher levels of hypoxia, indeed agreeing with the referenced experimental observations [1,16].

Cell-to-cell adhesion and its aggregative vs stabilising effect
Cell-to-cell adhesion, on the other hand, does not seem to play an important role in cluster formation or cluster topology. While small values of the cell-to-cell adhesion coefficient have an negligible effect on the simulation results overall, larger values show at least a transient impact on cluster formation and the cluster topology arising. While this weak effect may seem counter-intuitive considering previous single-cell vasculogenesis works [17,69,86,95], these models do not include both chemotaxis and degradation and thus do not predict cluster formation as an intermediate step in network formation. In view of these works, however, we cannot exclude that aspects of the final network topology, e.g. chord thickness, depend on the amount of cell-to-cell adhesion. Furthermore, according to Blatchley et al. [16], cell-to-cell adhesion is not related to cluster topology, which is coherent with our model predictions, but is responsible for cluster stability. The mathematical study of cluster-based vasculogenesis would therefore benefit from a more suited description of cell-to-cell adhesion, as we have seen in Section 3.3 that our modelling choice allows us to capture the aggregative effect of cell-to-cell adhesion, but not its role in cluster stabilisation. This could perhaps be achieved by, in addition to modelling cell-to-cell and cell-to-matrix adhesion nonlocally, modelling diffusion nonlinearly [23,72]. Alternatively, one could modify the volume exclusion term (i.e. definition (3)) to explicitly model cell-adhesion molecules, which would allow us to study the effect of cell-adhesion-mediated saturation of chemotaxis [68,70,103]. Otherwise, the stabilising effect of cell-to-cell adhesion could be more easily achieved by adopting an individual-based or hybrid modelling approach [5,103,113,114], as done in previous CP models of single-cell vasculogenesis [17,69,70,86,96,95,107].

Experimental research perspectives
The results reported in this paper suggest that in order to gain a comprehensive understanding of EPC cluster formation during the early-stages of cluster-based vasculogenesis, it is necessary to consider both ECM degradation and chemotaxis. So far, these processes have only been investigated separately in experimental studies [1,16] but it would be interesting to examine them in both normoxic and hypoxic conditions, in view of our in silico results (cf. Section 4.1). Model validation would also allow for a more accurate baseline parameter set. In fact, while the role played by ECM degradation and the size of clusters in our numerical simulations well matches the experimental observations reported by Blatchley and coworkers [16], this -and future -mathematical models would benefit from the estimation of parameter values by fitting the model to experimental data in order to obtain a better match of the timescale of 2D cluster formation. Note, for instance, that the ten-fold increase in α m which results in the observation of well-defined clusters at t = 20 (see supplementary Figure S7) might be justified by the experimental conditions reported in [16] (see A). In view of the results reported in Section 3.5 in absence of matrix remodelling, it would also be interesting to experimentally investigate the role ECM remodelling may plays in cluster formation and cluster size. This could lead to interesting observations regarding the empirical differences of cluster-based vasculogenesis during in vitro and in vivo studies, and suggest interesting focus points for therapeutic intervention. In order to achieve this, new in vitro experimental set ups might be required, in order to avoid the remodelling-resembling effects of cells falling towards the bottom of the 3D hydrogel during the cluster-formation process.

Mathematical research perspectives
In this work we have conducted a parametric analysis of the model in order to elucidate potential model behaviour. Whereas this already gave meaningful information to discuss the importance of the various processes involved in cluster-based vasculogenesis there are more advanced mathematical tools available for the analysis of the dependence of the model outcome on the model parameters. In particular, once accurate ranges of parameter values, possibly together with a distribution of the values over the range, have been estimated for various parameters from data spanning different levels of hypoxia, it would be of significant interest to conduct a global (or partially global) sensitivity analysis [65,85,87] of the cluster width and compactness to parameter variability. In fact, while the conducted parametric analysis highlighted the role played by each single parameter as it deviates from its baseline value, the global sensitivity analysis would provide a tool to investigate the effect of combined variations of multiple or all parameters over ranges indicative of different levels of hypoxia. This global approach is also significant given the nonlinear nature of the model [94,93]. Upon revision of the modelling strategy adopted to describe cell-to-cell adhesion, a global sensitivity analysis could be conducted in 2D in relation to the cluster compactness (C), area (W ) and elongation (see Section 2.5). Moreover, it would be particularly interesting to conduct a weakly nonlinear analysis [18,30,44,51] of the 2D model for a quantitative description of how chemotaxis is responsible for cluster topology. For instance, considering the results reported in Figure 10, there might be a threshold value of the chemotactic sensitivity below which we do not observe 2D clusters.

Modelling research perspectives
In addition to the modelling research perspectives discussed in Sections 4.1 and 4.2, interesting investigations could be conducted by considering the effect of persistence of motion or cell-matrix mechanical interactions.

The possible interplay between degradation and persistence of motion
Previous CP models of single-cell vasculogenesis predicted cluster formation in absence of cell elongation [67], due to chemotaxis. Despite their exclusion of ECM degradation and proliferation, we expect such result to capture the role of elongation during cluster-based vasculogenesis, given that EPC motion during the early stages of the process is mostly amoeboid-like, and we observe cell polarisation only in the later stages [16].
On the other hand, persistence of motion, indicated to prevent cell clusterisation by PEC models of single-cell vasculogenesis [111], was observed to be enhanced during amoeboid-like migration [98]. It would therefore be interesting to investigate whether clusters would still form in this modelling framework with the inclusion of persistence of motion. In particular note that Tosin et al. [111] included in their momentum equation a drag force generated by cells moving on the ECM and observed that lower cell-matrix adhesion resulted in cell clusterisation. We might thus expect ECM degradation to play a key role in lowering cell adhesiveness to the ECM in such modelling framework, thus leading to clusters.

Modelling the late-stage dynamics of cluster-based vasculogenesis
We plan to extend the theoretical investigation of cluster-based vasculogenesis to later stages of this process, during which EPCs increase their interaction with the ECM and bridge clusters, forming the vascular network. At such stage the mechanical interaction between the cells and ECM becomes non-trivial and we are therefore going to consider a mechanochemical model similar to those previously proposed to study the late stages of single-cell vasculogenesis [64,74,111]. Various works in the current literature already address the formation of sprouts from existing clusters [17,70,106,107] and prior single-cell vasculogenesis models investigating the formation of a vascular network, rather than clusters, may still be relevant to these late stages. For instance, while cell-to-cell adhesion may saturate chemotaxis [70], VEGF gradients may still be strongly perceived at the cluster boundaries. This is especially the case if ECM-bound VEGF was the main signalling molecule, and in this case we would indeed expect steep VEGF gradients at the cluster boundaries [67]. Perhaps such a scenario can explain the change in cell morphology. The role cell elongation plays in the formation of a final well-defined network has already been extensively addressed by CP models of single-cell vasculogenesis [17,69,67,81,86,95,116] Blatchley et al. [16] indeed reported an elongated cell morphology during sprouting, as well as an increased mechanical interaction with the ECM. The combination of cell traction on the ECM and strain-dependent movement of cells would also suffice in explaining bridging between clusters [3,64,74,76]. Nevertheless mechanochemical models, as well as CP models including cell traction [86,116], generally predict that higher ECM stiffness inhibits network formation, which is in contradiction with the experimental observations presented in [16]. Existing single-cell vasculogenesis models of late-stage dynamics, however, do not include ECM degradation, a key element for the study of cluster-based vasculogenesis. On the other hand, ECM degradation has been shown to have an important role during sprouting angiogenesis [17,33,50,95,112], a result which could very well translate into an active role in cluster-based vasculogenesis -although the interplay between ECM degradation and cell-matrix interactions may be particularly complex. Upon formulation of a mathematical model which accurately predicts cluster-based vascular network assembly, this could be used to investigate the determinants of network size and topology. For instance, one might explore whether the size of clusters or cell-to-cell adhesion [70,86] will affect tube diameters, or whether VEGF diffusion and decay rates determine chord length [4,3,38,98], and how these are affected by matrix stiffness.

A Parameter values
As already indicated in Section 2.2, we use L = 0.1 cm as characteristic length scale, in accordance with previous vasculogenesis works [64,98] and for easy visual comparison with the experimental results reported by Blatchley and coworkers [16]. We then take reference time scale τ := L 2 /D, where D is a characteristic diffusion coefficient D ∼ 10 −6 cm 2 s −1 [19], resulting in a reference time scale τ = 10 4 s.

Endothelial progenitor cells (or endothelial cells)
The reference cell density is chosen to be N := n M = ϑ −1 1 and we take ϑ 1 to be the average volume occupied by an EC. Rubin and coworkers [89] measured the average ECs volume during different phases of the cell cycle, registering values in the range of 800 − 1800 µm 3 , with a predominance of measurements around 1000 µm 3 . Hence we take ϑ 1 = 10 −9 cm 3 /cell, corresponding to an average cell diameter of approximately a = 10 −3 cm. Measures and estimates of the diffusion coefficient of ECs fall in the range 10 −6 -10 −12 cm 2 s −1 [3], so we take D n = 10 −9 cm 2 s −1 . We consider the chemotactic sensitivity coefficient estimated by Jain and Jackson [52], corresponding to χ = 1.4 × 10 −7 cm 5 ng −1 s −1 . Sen et al. [97] estimated a maximum cell surface radius, upon morphological changes to better adhere to the underlying gel, of about 50 µm. Given cell-to-cell and cell-to-matrix adhesion occurs via adhesion molecules on the cell surface, we take the sensing radius R = 0.5 × 10 −2 cm. While the nonlocal term introduced in equations (5) and (6) allows us to well consider cell-to-cell and cell-to-matrix adhesion dynamics at tissue level, these are the result of smaller scale dynamics between adhesion molecules and receptors on cell surfaces -vid. for instance [2,13,39] -of which our modelling choice is a simplification. As a result, good estimates for cell-to-cell and cell-to-matrix coefficients S nn and S nρ are currently lacking. We therefore consider nondimensional values chosen in [42] for our baseline parameter set which correspond to S nn = 10 −16 cm 5 s −1 and S nρ = 10 −6 cm 2 nM −1 s −1 respectively, acknowledging that model fitting to experimental data is required. Kinev and coworkers [55] reported dubling times of non-irradiated ECFCs -the same class of EPCs employed by Blatchley et al. [16] -of about 19.5 hours, estimated from measured growth rates assuming exponential growth. Following their calculations, this corresponds to proliferation rates of about p = 10 −5 s −1 .
Extracellular matrix We use a reference matrix density of P = 10 −1 nM [5,6,110] and define the parameter ϑ 2 := P −1 . We take the ECM degradation rate by MMPs proposed by Kim and Friedman [54], i.e. γ = 9 × 10 5 cm 3 g −1 s −1 . ECM remodelling is a complex process involving a variety of cells and molecules [24,32,61,105], so the remodelling term introduced in equation (10) is an oversimplification of the underlying dynamics. Therefore the lack of experimental values or estimates for the remodelling rate µ is not surprising, and we take the nondimensional value 0.2 as similarly considered in [34,37,42] for our baseline parameter set. Under the defined nondimensional parameters this corresponds to the dimensional value µ = 0.2 × 10 −5 nM s −1 .
Matrix degrading enzyme (MMP) Blatchley et al. [16] reported concentrations of MMP-1 in the range 1 − 100 µg ml −1 , so we take the intermediate concentration as reference MMP density, i.e. M = 10 µg cm −3 . We let the diffusion coefficient for the MMP be given by D m = 8 × 10 −9 cm 2 s −1 , which was experimentally determined by Saffarian and coworkers [92], although diffusion rates have been observed in the range 10 −10 −10 −8 cm 2 s −1 [27,59]. While MMPs secretion rates by cells have been reported in a variety of works, these fail to provide parameter values in the appropriate unit of α m in this model. For instance in [59] secretion rates considered varying between 0.01 − 0.5 s −1 but the associated MMP concentration is unspecified, while in [90] the estimated rate is of 5.75 × 10 −10 and 4.44 × 10 −10 g cm −3 secreted by stromal cells and macrophages respectively, without an indication of the considered time frame. We therefore here consider that molecular dynamics are generally faster than cellular ones and therefore the observed average MMP concentration satisfies (11) at equilibrium, i.e. m = α m n/λ m . Under the chosen reference values N , M and λ m , this corresponds to an MMP production rate of α m = 0.5 × 10 −12 µg s −1 per cell. Note that the resulting nondimensional parameter value is close to those chosen in previous mathematical models -vid. for instance [6,34,37,42]. Nonetheless, we will also consider higher values of α m , since Deem and Cook-Mills [35] have reported up to a 4-fold increase in MMP production in the presence of reactive oxygen species, and in addition MMP production levels have been shown to be significantly upregulated in human cancers [102]. Finally, Kim and Friedman [54] estimated -from its half life -the decay rate of MMP λ m = 5 × 10 −5 s −1 .

Chemoattractant (VEGF)
We take the reference VEGF density to be C = 20 ng cm −3 , in the range of values generally considered in in vitro set ups -vid. for instance [45,60,98]. Measures and estimates of the diffusion coefficient of the chemoattractant, usually identified as VEGF, are in the range 10 −6 − 10 −9 cm 2 s −1 [3,38,70,71,98,103], so we take D c = 10 −7 cm 2 s −1 . Yen and coworkers [118] reported a VEGF secretion rate of 0.068 molecules cell −1 s −1 , which in combination with Avogadro's number (6.022 × 10 23 molecules per mole) and a molecular weight of 45kDa [118] results in a VEGF production rate of about α c = 0.5 × 10 −11 ng s −1 per cell. Serini et al. [98] reported the half-life of VEGF-A to be approximately 64 minutes, corresponding to a decay rate of about λ c = 2.7 × 10 −4 s −1 , in line with the values chosen in [70,103].
A nonlocal partial differential equation model of endothelial progenitor cell cluster formation during the early stages of vasculogenesis * Supplementary Information In this document (supplementary information) we report supplementary material to the contents of the above named publication (main publication). In particular we include detailed linear stability analysis calculations in Section S1, provide supplementary figures in Section S2, and explain a detail related to the numerical treatment of the zero-flux boundary condition of the nonlocal term, in Section S3 of this document.

S1 Linear stability analysis
We here conduct a linear stability analysis (LSA) on the model equations. To understand the role played by matrix degradation on the evolution of the solution of our model -and therefore its role in the early stages of cluster based vasculogenesis -particular emphasis will be put on the differences between the case in which γ > 0 and γ = 0. First we look at the full nondimensional system (15) with saturation effects, that is the one making use of definition (16) for the function f (n, ρ). We will also consider the simplified version of the model under definition (17) instead, corresponding to the absence of saturation effects, for a better understanding of the results of the LSA. The analysis is carried out for the 1D problem for simplicity, but in section S1.5 we will briefly consider the corresponding 2D problem.

S1.1 Spatially homogeneous steady states
The spatially homogeneous steady states of the system (15)v = (n,ρ,m,c) satisfy giving eitherv = (0, 1, 0, 0) orv = (1, 0, α m /λ m , α c /λ c ) . Each spatially homogeneous steady state is therefore in the formv with eithern = 0 orn = 1. The first one corresponds to the absence of cells, with the whole volume occupied by the ECM. On the other hand, the second one corresponds to the case in which the cells have completely degraded the ECM and are occupying the whole volume. Note that in absence of matrix degradation, i.e. for γ = 0, the spatially homogeneous steady states are still in the form (S2), this time with 0 ≤n ≤ 1. We therefore have, in addition to the two states already described, infinitely many steady states in which both cells and ECM are present, filling up the volume. These additional steady states might be referred throughout this section as "intermediate" steady states, as they correspond to intermediate values ofn, namely 0 <n < 1.

Remark 1:
We are going to assume that n + ρ ≤ 1 throughout the rest of the section. That is indeed satisfied by the steady states and we justify it through the biological argument that random perturbations to these states are likely to arise naturally if there is space available, while changes that would imply n+ρ > 1 are much less likely as they would require high energy expenses to oppose a high pressure environment.
Mathematically this allows us to avoid complications that arise when introducing perturbations in the typical ansatz (vid. below) due to definitions (3) and (10). In fact, under this assumption the term (1 − n − ρ) + = (1 − n − ρ) ≥ 0 for all perturbations allowed. Remark 2: During the linear stability analysis reported in the following sections, we are going to use the fact thatn +ρ = 1 and γρn = γρm = γρc = 0 for all spatially homogeneous steady states, as concluded above, to further simplify calculations. Remark 3: In the following we are going to consider biologically significant steady states, that is we are going to only investigate the stability of the steady states with n > 0. The stability of the cell-free steady state n = 0 could be investigated by considering the dynamics along a center manifold about the steady state, as linear stability analysis calculations reported below are inconclusive for this case. Such investigation, however, goes beyond the scope of this study and we leave the mathematical details to the interested reader, referring them to [?, ?]. However, simple numerical simulations indicate that this cell-free steady state is unstable.

S1.2 Stability under spatially homogeneous perturbations
Introducing a small spatially homogeneous perturbation v =v+ṽ(t), with |ṽ| 1, in (15) and linearising leads to the following system for the perturbationṽ(t): 2 Assuming small perturbations in the formñ,ρ,m,c ∝ exp (σt), the system (S3) can be rewritten as For a non-trivial solution we require det M = 0, leading to the dispersion relation which can only be satisfied for Re(σ) < 0, indicating the biologically significant spatially homogeneous steady state is stable under spatially homogeneous perturbations. On the other hand, in the absence of matrix degradation, i.e. for γ = 0, we obtain the corresponding version of the dispersion relation (S5), which can also be satisfied for σ = 0, thus linear stability analysis is inconclusive. In this case, under the assumption specified in Remark 1, the additional spatially homogeneous steady states are actually unstable to spatially homogeneous perturbations, as can be quickly verified numerically, and the system will converge to a nearby steady state depending on the initial conditions of n and ρ, as well as the values of parameters p and µ.

S1.3 Stability under spatially inhomogeneous perturbations
Let us now introduce spatially inhomogeneous perturbations v =v +ṽ(t, x), with |ṽ| 1. After linearisation system (15) gives We immediately notice that, as all steady states satisfy volume filling conditions, chemotaxis in presence of saturating effects does not play a role in the dynamics of small perturbations from any of these states.
Assuming small perturbations in the formñ,ρ,m,c ∝ exp (σt + ikx), note first of all that (S7) can be rewritten as We will make use of the notation for which we have that w 1 (k) ≥ 0 for all k ∈ R. Then, for small perturbations in the form introduced above, the system (S6) can be rewritten as with w 1 (k) defined in (S9). This immediately implies A(k) ≥ 0 for all steady states, those with and without matrix degradation. The dispersion relation σ(k 2 ) can then be obtained by imposing det M = 0 for a non-trivial solution to (S10), and in order for the steady states to be unstable to spatially inhomogeneous perturbations -and patterns to arise -we require Re σ(k 2 ) > 0 for some k 2 ∈ R. This, in its more general form, is given by It is clear from (S12) that Re σ(k 2 ) < 0 for all k 2 . Therefore all steady states considered are also stable to spatially inhomogeneous perturbations. We claim this to be due to saturation effects introduced in (3), while we may still expect chemotaxis and cell-to-cell or cell-to-matrix adhesion to play an important role for cell aggregation when the initially perturbed states are far from being volume filling. Let us explore the matter further in the next section.

S1.4 Considerations in absence of saturation effects
We here consider the pattern formation potential of the model in absence of saturation effects, that is we make use of definition (17) for the function f (n, ρ), which influences both the chemotactic sensitivity of cells and their adhesion velocity. All results and considerations made in sections S1.1 and S1.2 on spatially homogeneous steady states and their stability under spatially homogeneous perturbations still hold.
Introducing spatially inhomogeneous perturbations as in section S1.3 leads to the following linearised system alternative to (S6): in which, for the 1D problem, A[v +ṽ(t, ·)] is given by x + rη(j)) + S nρρ (t, x + rη(j)) dr (S14) with Γ(r) given by the corresponding definition in (8). Notice that under these modelling assumptionsnamely no saturation effects included in the adhesion velocity -we do not need to assume n + ρ ≤ 1 for analytical tractability, and therefore have no crafty restrictions on the perturbations introduced. On the contrary, the contents of Remark 2 and Remark 3 still hold. Then, assuming small perturbations in the formñ,ρ,m,c ∝ exp (σt + ikx), (S14) can be rewritten as with w 1 (k) ≥ 0 defined in (S9). We can then rewrite (S13) as Mṽ = 0, with M given by where We have both A n (k) ≥ 0 and A ρ (k) ≥ 0 for all k ∈ R. Comparing these to (S10) and (S11), we see that in absence of saturation effects the contributions from cell-to-cell and cell-to-matrix adhesion are decoupled and both give a negative contribution to M. Furthermore, the contribution from chemotaxis does not vanish for any of the biologically significant steady states (i.e.n > 0). From this system the dispersion relation σ(k 2 ) satisfies (σ + D m k 2 + λ m ) (σ + D n k 2 ) σ + γm + µ + pn − A n (k) σ + γm We now have that any biologically significant steady state may be unstable to spatially inhomogeneous perturbations for strong enough cell-to-cell adhesion and/or chemotaxis, i.e. for large enough S nn and α c and/or χ. In addition, the magnitude of cell-to-cell adhesion A n (k) in (S17) will be larger for smaller values of the sensing radius R and larger values ofn, and the contribution from chemotactic movement will also increase for larger values ofn. Therefore in absence of saturation effects steady states corresponding ton > 0 may be unstable to spatially inhomogeneous perturbations and patterns may form. This suggests that perturbed initial conditions far from being volume filling, so that saturation effects do not play a big role, might result in cell aggregation thanks to cell-to-cell adhesion and chemotaxis, as long as the initial cell density -interpreted as a small perturbation from a spatially homogeneous steady state -is large enough for these dynamics to play a significant role.

S1.5 Extension of LSA results to the 2D problem
We here consider how the results obtained so far change in the 2D problem. Conclusions drawn in sections S1.1 and S1.2 remain unchanged.
Considerations in presence of saturation effects. When introducing spatially inhomogeneous perturbations v =v +ṽ(t, x), after linearisation A[v +ṽ(t, ·)](x) for the 2D problem is given by x + rη(θ)) +ρ(t, x + rη(θ)) dθ dr (S19) with Γ(r) given by the corresponding definition in (8). Assuming small perturbations in the form n,ρ,m,c ∝ exp (σt + ik · x), with k = (k 1 , k 2 ) ∈ R 2 , then (S19) can be rewritten as where we have defined, with notation that will become more intuitive in a few steps, We eventually obtain the system (S10), in which k 2 = |k| 2 = k 2 1 + k 2 2 and, instead of A(k) as in (S11), we have A(k) given by Figure S1: Plot of the quantity −k · Θ R (k), where k = (k 1 , k 2 ) and Θ R (k) is defined in (S21) for nondimensional R = 0.05. We clearly have −k · Θ R (k) ≥ 0 for all k ∈ R 2 , ensuring the real part of A(k) in (S23), as well as those of A n (k) and A ρ (k) in (S26) are non-negative. The integral in (S21) was evaluated numerically using the trapezoidal rule for numerical integration.
with Θ R (k) and Θ I (k) defined in (S21) and (S22). Despite the complexity of these integrals, we evaluated numerically the quantity −k · Θ R (k) -see Figure S1 -which allowed us to conclude that, just like in the 1D case, Re A(k) ≥ 0 for all k ∈ R 2 . For this problem we therefore obtain the corresponding version of equation (S12) for the dispersion relation and reach the same conclusions as in section S1.3 for the full system including saturation effects.
Then assuming small perturbations in the formñ,ρ,m,c ∝ exp (σt + ik · x), with k = (k 1 , k 2 ) ∈ R 2 , (S24) can be rewritten as with Θ R (k) and Θ I (k) defined in (S21) and (S22). We eventually end up with a system Mṽ = 0, with M defined as in (S16), in which again we have k 2 = |k| 2 = k 2 1 + k 2 2 and, instead of A n (k) and A ρ (k) as in (S17), A n (k) and A ρ (k) defined by Again, thanks to the numerically evaluated quantity plot in Figure S1, we have that Re A n (k) ≥ 0 and Re A ρ (k) ≥ 0 for all k ∈ R 2 , and can therefore reach analogous conclusions to those drawn in section S1.4 for the corresponding 2D problem. Figure S2 contains the plots up to t = 400 of the same solutions plot in Figure 6 in the main publication, in order to show long-time asymptotics of such solutions. Figure S3 contains plot analogous to those in Figures 7 and 8 in the main publication, with additional information complementing the results presented in Section 3.4 of the main publication. In particular the plots display how choosing parameter values of different orders of magnitude compared to those chosen for the baseline parameter set (see Table 1 in the main publication) will influence average cluster width W and compactness C, for the parameter D n (Figure S3a), S nn ( Figure S3b), γ ( Figure S3c), D m ( Figure S3d) and D c ( Figure S3e). Figure S4 contains plots of the two-dimensional simulations under the baseline parameter set shown in the second row of Figures 9, 10 and 11 of the main publication, displaying the solution at additional times, showing the long-term behaviour of the solution. Figure S5 contains plots of the two-dimensional simulations under the baseline parameter set in absence of matrix degradation (i.e. γ = 0), starting from the same initial conditions (cf. first panel in first row of Figure S4). Figure S6 contains plots of the two-dimensional simulations under the baseline parameter set (second row) and under different values of the MMP secretion rate α m , starting from the same initial conditions (cf. first panel in first row of Figure S4). Figure S7 contains plots of the two-dimensional simulations under the baseline parameter set (second row) and under different values of the VEGF secretion rate α c , starting from the same initial conditions (cf. first panel in first row of Figure S4). Figure S2: First row: Plots of the cell density n(t, x) up to t = 400 obtained solving the system (15), together with definitions (6), (8) and (16), initial conditions (13) and (14) 1 , complemented with zero-flux boundary conditions, in absence of cell-to-cell adhesion, i.e for S nn = 0: under the baseline parameter set (first column), in absence of matrix degradation, i.e. for γ = 0 (second column), in absence of chemotaxis, i.e. χ = 0 (third column), and in absence of both chemotaxis and matrix degradation, i.e. χ = γ = 0 (fourth column). Second, third and fourth rows: Same as first row but in the presence of cell-to-cell adhesion, with S nn = 0.1 (second row), S nn = 1 (third row) and S nn = 10 (fourth row) respectively.  (20) and (21), under variations of certain parameters from the baseline parameter set (BPS), in Table 1 of the main publication. In (a) we have boxplots of W and C measured on the numerical solution of the system (15) at t = 50, for D n taking its value in the BPS (center), one order of magnitude higher (left) and one order of magnitude lower (right) than the one in the BPS. Each boxplot collects data from 100 simulations under randomised initial conditions (13) and (14) 1 . In (b)-(e) we have the same as in (a) but varying parameters S nn (b), γ (c), D m (d) and D c (e). Figure S4: Plots of the cell density n(t, x) obtained solving the system (15), together with definitions (6), (8) and (16), initial conditions (13) and (14) 2 , complemented with zero-flux boundary conditions, under the parameter choices reported in Table 1. The solution is plotted at time t = 0, 20, 40, 60, 80 (first row, left to right), t = 100, 120, 140, 160, 180 (second row, left to right) and t = 200, 220, 240, 260, 500 (third row, left to right). Figure S5: Plots of the cell density n(t, x) obtained solving the system (15), together with definitions (6), (8) and (16), initial conditions (13) and (14) 2 , complemented with zero-flux boundary conditions, under the parameter choices reported in Table 1, except for γ = 0. The solution is plotted at time t = 20 (first panel), t = 40 (second panel), t = 80 (third panel) and t = 140 (fourth panel). Figure S6: First row: Plots of the cell density n(t, x) obtained solving the system (15), together with definitions (6), (8) and (16), initial conditions (13) and (14) 2 , complemented with zero-flux boundary conditions, under the parameter choices reported in Table 1, except for α m = 0.05. The solution is plotted at time t = 20 (first panel), t = 40 (second panel), t = 80 (third panel) and t = 140 (fourth panel). Second and third row: Same as first row, except for α m = 0.5 (second row) and α m = 5 (third row). Figure S7: First row: Plots of the cell density n(t, x) obtained solving the system (15), together with definitions (6), (8) and (16), initial conditions (13) and (14) 2 , complemented with zero-flux boundary conditions, under the parameter choices reported in Table 1 S3 Toeplitz-to-circulant embedding in the evaluation of the nonlocal term with zero-flux boundary conditions Consider the nonlocal term as given in Equation (5) in Section 2.1.1 of the main document, that is

S2 Supplementary figures
together with our realisation of zero-flux boundary conditions, that is g v(t, x) = 0 if x is outside of the spatial domain Ω. A detailed account of the accurate and efficient numerical treatment of (S27) with periodic boundary conditions, also for the corresponding spatially two-dimensional case, is given in [?]. In the case of periodic boundary conditions, an efficient numerical scheme is possible if the spatial domain is covered by a uniform grid with N grid cells because then, on each grid cell interface x, the same discretization of the nonlocal term in x can be applied. This gives rise to an evaluation of the nonlocal term on all N grid cell interfaces, recall the periodic boundary condition, simultaneously. We denote the vector of these N G for our example have the form We can now embed T into a circulant matrixC ∈ R N +2,N +2 of minimal size in the following waỹ Note that the actual dimension ofC depends on the upper and lower bandwidth of matrix T . Finally, definingÃ we can computeÃ =CG using FFT techniques as in the periodic boundary condition case and recover the desired vector A as the first N − 1 entries ofÃ. Thanks to the Toeplitz-to-circulant embedding described here, also nonlocal models with zero-flux boundary conditions can be simulated efficiently. The techniques outlined here carry over to the spatially twodimensional setting, where a block-Toeplitz matrix with Toeplitz blocks is embedded in a block-circulant matrix with circulant blocks and 2D-FFT is applied, see [?].