# Modelling of chemotactic sprouting endothelial cells through an extracellular matrix

^{1}Department of Condensed Matter Physics, University of Barcelona (UB), Barcelona, Spain^{2}Electronics and Biomedical Engineering, University of Barcelona (UB), Barcelona, Spain^{3}Biomaterials for Regenerative Therapies, Institute for Bioengineering of Catalonia (IBEC), The Barcelona Institute of Science and Technology (BIST), Spain^{4}Institute for Research in Applied Mathematics and Systems, National Autonomous University of Mexico , Mexico City, Mexico^{5}Institute of Physics, National Autonomous University of Mexico, Mexico City, Mexico^{6}Institute of Nanoscience and Nanotechnology (IN2UB), University of Barcelona (UB), Barcelona, Spain

Sprouting angiogenesis is a core biological process critical to vascular development. Its accurate simulation, relevant to multiple facets of human health, is of broad, interdisciplinary appeal. This study presents an *in-silico* model replicating a microfluidic assay where endothelial cells sprout into a biomimetic extracellular matrix, specifically, a large-pore, low-concentration fibrin-based porous hydrogel, influenced by chemotactic factors. We introduce a novel approach by incorporating the extracellular matrix and chemotactic factor effects into a unified term using a single parameter, primarily focusing on modelling sprouting dynamics and morphology. This continuous model naturally describes chemotactic-induced sprouting with no need for additional rules. In addition, we extended our base model to account for matrix sensing and degradation, crucial aspects of angiogenesis. We validate our model via a hybrid *in-silico* experimental method, comparing the model predictions with experimental results derived from the microfluidic setup. Our results underscore the intricate relationship between the extracellular matrix structure and angiogenic sprouting, proposing a promising method for predicting the influence of the extracellular matrix on angiogenesis.

## Introduction

The ability of cells to migrate is essential for fundamental biological processes involving vasculature development, angiogenesis and wound healing Qiu and Hirschi (2019); Friedl and Gilmour (2009). Specifically, migrating endothelial cells (ECs) are crucial in the formation of a functioning vascular network and angiogenesis. EC constitutes the inner cellular lining of capillaries, arteries and veins. ECs need to adapt, forming new vessels responding to tissue injury or lack of nutrient conditions. To produce a functional vasculature, growth and transcription factors interplay to balance angiogenesis and EC migration Jain (2003); Krüger-Genge et al. (2019); Carmeliet and Jain (2011).

Sprouting angiogenesis, or the penetration of vessels into avascular tissue, is critically dependent on the interactions between ECs, pericytes, and the extracellular matrix (ECM) Marchand et al. (2019); Blocki et al. (2018). The ECM plays an integral role in supporting and guiding the cellular sprouting process. Initially, endothelial cells extend filopodia to probe the surrounding ECM environment, guided by gradients of bioactive molecules De Smet et al. (2009). The ECM, in turn, provides a physical scaffold for these exploratory cells and presents a rich source of sequestered growth factors that can be liberated by enzymatic degradation Schultz and Wysocki (2009). Key among the molecular regulators of angiogenic sprouting is the vascular endothelial growth factor (VEGF), which promotes the proliferation and migration of endothelial cells. In addition to VEGF, the Delta-like ligand 4 (Dll4)-Notch signalling pathway plays a critical role in regulating the selection and behaviour of the sprouting cells Pitulescu et al. (2017). The dynamic interplay between cells and the ECM, conditioned by the liberation of sequestered proteins and bioactive molecules, lies at the heart of angiogenic sprouting. This complex relationship not only shapes cellular behaviours during angiogenesis but also contributes to the stabilisation and maturation of the emerging vascular network. Furthermore, integrins, a family of cell surface receptors, are central to cell-ECM interactions, mediating cell adhesion to the ECM and transducing biochemical and biomechanical signals from the ECM into the cell Kanchanawong and Calderwood (2023). Recognizing the pivotal role of integrins and bioactive molecules in mediating the dynamic interplay between cells and the ECM, the significance of the ECM extends beyond a mere structural scaffold. Indeed, it serves as an active participant in cellular processes, functioning as a complex and adaptable microenvironment. Through the ECM, cells can migrate, proliferate, sense, communicate with other cells, and even remodel this same ECM in a sophisticated and spatio-temporal synchronisation Lutolf and Hubbell (2005). The reciprocal cell-matrix interactions and the different species diffusing along the ECM determine many of the cell behaviour. Because of the strong interaction between cells and their surrounding, the fabrication of artificial ECM has become of tremendous interest and challenge.

Hydrogels, in many forms, are widely used as extracellular matrix materials for *in vitro* experiments Vila et al. (2020); Krishnamoorthy et al. (2019); Blache and Ehrbar (2018); the use of hydrogels is so far the only option, given their biocompatibility, softness and water content Caló and Khutoryanskiy (2015). Nowadays, tissue engineering can promote angiogenic growth factor release and vascularisation to irrigate tissues. However, still, there is a lack of control and synchronisation between fabricated ECMs, vasculature structures and cell fates. Therefore, theoretical studies and computational modelling add value to understanding and elucidating how the structure of hydrogels, or artificial extracellular matrices, can influence cell dynamics Vats and Benoit (2013); Ma et al. (2021); Zhang K. et al. (2021); Elosegui-Artola (2021). Fibrin-based hydrogels have been used as ECM analogues thanks to their feasibility of tuning their mechanical properties by simply changing concentrations or gelation conditions.

In this study, we have performed experiments to study the advancement of ECs through a large-pore, low-concentration fibrin-based hydrogel matrix due to a chemotactic factor gradient concentration. This hydrogel matrix generates a porous landscape for cells to migrate, replicating their natural environment, i.e., acting as an ECM. The mechanical properties of the hydrogel are characterised at the macro-scale in the Supplementary Material, *“Characterisation of the hydrogel structure”*. In (Figure 1A) we show the experimental design. The microfluidic platform consists of a central chamber for the hydrogel and two side channels, with two reservoirs each. Reservoirs are filled up with different cell media volumes to generate flow inside the side channels due to the difference in hydrostatic pressure. At one side of the central chamber, we have grown an endothelium using rat endothelial progenitor cells (rEPCs). On the other side-channel, we placed a chemotactic factor, VEGF, that diffuses through the ECM. We see in (Figures 1B, C) the initial distribution of cells and their response to the chemotactic factor gradient after 48 h.

**FIGURE 1**. Schematic representation of the experimental setup and results of cell migration. **(A)** Scheme of the microfluidic platform design for angiogenic cell migration within the central chamber (1300 *μm* wide, 8800 *μm* long and 150 *μm* high) filled up with fibrin-based hydrogel (violet). Two lateral channels are connected to the central chamber through a trapezoidal structure (dark grey) with deposits at their ends. One side channel contains rat endothelial progenitor cells (rEPC) (grey), and the opposite side channel is connected to reservoirs with angiogenic factor concentration (vascular endothelial growth factor, VEGF) (yellow). **(B)** Fluorescence microscopy image of the cells inside the microfluidic system without angiogenic factor after 48 h. **(C)** Fluorescence microscopy image of the cells after 48 h under the influence of the angiogenic factor gradient. Endothelial cells migrate from an initial endothelium or isle, at one side of the channels, towards the higher concentration of VEGF, forming a finger-like structure. The scale bar in **(B,C)** represent 50 *μm*. Cells were stained for F-actin (red), and cell nuclei (blue). Images **(B,C)** depict the trapezoidal structure marked with white dashed lines, rotated 90° from image **(A)**.

ECs tend to reach a quiescent mode forming an endothelium. This endothelium is formed at the border side of the hydrogel, mimicking the ECM in the side channel. We see a picture of the artificial ECM in (Figure 2A). Cells remodel and advance through the pores of the ECM, physically enlarging the pores and chemically degrading the extracellular matrix. Because of the changes cells produce in the ECM matrix liberating proteins, ECM-cell interactions become more intricate and relevant. The ECM conditions the behaviour of cells and acts as a scaffold to stabilize the structure of the forming vasculature network.

**FIGURE 2**. Cell environment and modelled cell environment: confocal microscopy image of a section of the 3-D hydrogel using a Z-stack. **(A)** Confocal microscopy image of a region of 775* μm*^{2} of fibrin-based hydrogel used for cell migration; the image has been generated by NHS-rhodamine covalently linked to amine groups included in fibrin. *In-vitro* experiments use this hydrogel to study endothelial cell migration and angiogenesis. The colour scale represents the hydrogel concentration present at each pair of coordinates. The scale bar is 100 *μm*. **(B)** Bidimensional reproduction of the experimental cell environment, the hydrogel and chemotactic factor. This description is used as a base for our simulations. The *C*_{0} environment is defined as the sum of a linear gradient of a chemotactic factor distribution (*ν*_{c}), increasing from green to yellow, and the extracellular matrix profile (*ν*_{ECM}), in purple.

Angiogenesis has been shown to rely on a migratory “tip”-cell. A tip-cell guides migration and is followed by proliferative stalk-cells that elongate the sprout Carmeliet et al. (2009); Jakobsson et al. (2010). The most sensitive EC to the chemoattractant stimuli becomes a tip cell that guides the vessel sprout at the forefront. This tip cell forms numerous filopodia that scan the environment for angiogenic cues De Smet et al. (2009); Mayor and Etienne-Manneville (2016). We also see this collective behaviour in our experiments (Figure 1C). In this image we observe a collective of cells advancing a central cytoplasmatic finger shape that should become a capillary in time. The growth and dynamic evolution of this angiogenic sprouting is what we model and characterize in this study.

To study the behaviour of cells, as a collection, seen in (Figure 1C) we consider both the ECM structure and the chemotactic factor gradient present in our experiment. We model the chemotactic collective migrating phenomena conditioned by the ECM structure. First, we neglect the remodelling of the ECM through mechanical and chemical interactions as the investigation focuses on the response of cells to the presence of chemical gradients and ECM distribution; we assume the cell environment to be spatial and temporally static. Then, we proceed with a model extension to allow cells to remodel the cell environment by sensing it and degrading the ECM. Various computational approaches have been applied to cell migration, and angiogenesis Alert and Trepat (2020); Heck et al. (2015), but cell motility is an outstanding example of a moving boundary problem. For this reason, some authors have proposed a phase-field theory to model angiogenesis Travasso et al. (2011); Moure and Gomez (2021); Lima et al. (2014); Vilanova et al. (2017); Xu et al. (2017); Giverso and Ciarletta (2016); Santos-Oliveira et al. (2015). Phase-field methodology solves most of the challenges related to the numerical solution of moving boundary problems and is one current methodology that has shown promising results in modelling single and collective cell behaviour Santra et al. (2020). The dynamics of cells in the phase-field methodology are governed by a partial differential equation that minimises an energy functional at a low computational expense.

There are two main ways to model cells: continuum models, which treat cells as concentrations of cell types, and discrete models, such as agent-based or particle models, which consider individual cells.

Cellular Automaton (CAs), Cellular Potts Models (CPMs) and Phase-Field Models are all lattice-based models representing the system as a set of cells or a continuum. On the other hand, Agent-Based Models (ABMs) are not necessarily lattice-based and can be either discrete or continuum. ABMs and CAs both use a set of rules to update the state of the system. In CAs, the state of each cell is updated depending on the neighbouring cells; meanwhile, the behaviour of ABMs is determined by the collective interaction of each individual agent. On the other hand, CPM and Phase-Field models are driven by energy minimization. In CPM, cells are treated as discrete entities, a collection of mobile lattice sites; meanwhile, Phase Field Models represent a continuous field that gradually changes across a transition region. Phase Field models represent the interface between two materials or regions by a transition zone where the properties of the materials gradually change.

The agent-based and discrete models described models above, provide a good intuition of the processes involved in the cellular tip extension and migration but need to explore the influence of the surrounding environment, which scale is considerably larger. The main problem in using these discrete rules to capture the response of global systems with variable conditions and high interactivity is that they rely on something other than physical principles but statistical or heuristic rules that control the extension and branching phenomena. These limitations of discrete models regarding the underlying development mechanics are addressed with continuum modelling by using partial differential equations (PDEs) governing more global descriptions of the presence of endothelial cells. Therefore, in recent studies hybrid models, combining both approaches, have gained traction as it allow us to gain a more comprehensive understanding of the system. However, appropriate linking and calibration of such models should be performed, and parameters are only sometimes measurable.

CPM has been extensively used to study cell migration and proliferation in angiogenesis Graner and Glazier (1992); Szabó and Merks (2013) and has also been extended as a hybrid model to account for different cases such as cell heterogeneity and its impact on tumours Leschiera et al. (2022); Zhang et al. (2023) and mechanical cell-matrix interactions van Oers et al. (2014). Other hybrid ABMs and CA models have been proposed to study the interaction between tumour and immune cells Ghaffarizadeh et al. (2018); Norton et al. (2019); Hendrata and Sudiono (2019). A different approach, closer to the approach we have applied to this study, is the example of the model presented in Travasso et al. (2011), which is a hybrid model that combines an agent-based rule to simulate the behaviour of individual cells and a PDE to account for the continuous model to capture the evolution of capillaries and angiogenic factors. This same model has been extended and applied to numerous situations since it allows tailoring the vascular network and proliferation rate by adjusting its agent-based parameters. Vilanova et al. (2013) proposed a high-order phase-field approach using isogeometric analysis. This model was subsequently refined and extended by Xu et al. (2016, 2017); Vilanova et al. (2017); Xu et al. (2020); Moure and Gomez (2021). These studies demonstrate the potential of phase-field modelling to capture the complexity of intermediate-scale growth dynamics.

There is a gap when looking for formulations that combine the continuous description of intermediate-scale physics with continuous integration of the micro-environment components. This paper proposes a new approach integrating environmental chemical interactions with cells using the phase-field modelling spirit. Therefore, it can capture the dynamics at the represented level without using explicit or discrete strategies.

Here we elucidate the influence of the ECM based on phenomenological phase-field modelling. Through quantitative matching of experiment and theory, we first reveal the occurrence of cell migration driven by the distribution of chemotactic factors and the structural configuration of the ECM. In a self-consistent manner, such revelation reflects the coupling strengths of the constituent materials of the ECM and biological factors, allowing for their quantification in the model. As opposed to other models, we write a single equation to describe cell migration and couple the cell environment contribution into a single parameter, making it more straightforward and interpretable. Our study thus demonstrates an approach to developing a comprehensive and predictive understanding of physical coupling phenomena for a potentially broad family of ECM and the distribution of cells. In this model, we focus on the dynamics and evolution of the conformation of the experimental results shown in Figure 1C, considering the chemotactic factor and the hydrogel structure. We aimed to recreate the process of sprouting angiogenesis, which specifically involves the migration of specialised ECs. Our phase-field methodology, describing collections of cells, allows us to represent the cellular response to its environment, i.e., the ECM and the distribution of the chemotactic factors.

### Experiments

We perform the experiments using the angiogenesis platform design represented in (Figure 1A), which has been inspired by the microfluidic platform used in previously published experiments López-Canosa et al. (2021, 2022). Deposits on each side of the channel have different concentrations of chemotactic factors, generating a gradient across the ECM. We use the invasion of ECs into the central chamber, filled with a large-pore, low-concentration fibrin-based hydrogel with a mean pore size of 40*μm*, to study chemotactic cell migration. The microfluidic chip is designed to separate the different phenomena to isolate and investigate the behaviour of endothelial cells in response to chemical signals.

Details of the setup were published previously López-Canosa et al. (2022) and are further elaborated in the *“Microfluidic assay and cell culture preparation”* section of the Supplementary Material. Briefly, the hydrogel (fibrinogen and bovine thrombin mixture) is injected into the central chamber and left for cross-linking. Rat endothelial progenitor cells (rEPC) are introduced into one of the side channels of the chip. Then the chip is incubated for 45 min in an upright position, laying on the opposite side of the cells, to allow cell attachment to the central chamber. Finally, the reservoirs were filled with medium with 120 *μL* upstream and 90 *μL* downstream.

In these experiments, we use VEGF as a chemotactic factor, not only because it is a very well-known chemoattractant but also because angiogenesis is usually promoted by the expression and generation of a VEGF gradient that activates ECs Shibuya (2001). In our experiments, the VEGF gradient is generated solely by the diffusion of the molecules through the extracellular matrix López-Canosa et al. (2022).

After culturing the cells for 24 h, we introduce the chemotactic factor into the system to induce cell migration through the ECM. The reservoirs are refilled with medium and medium + VEGF every 24 h to keep the hydrostatic pressure difference between the reservoirs, which generates the flow in the side channels and allows us to have a constant gradient distribution as evidenced in López-Canosa et al. (2022), we see that the gradient of a 40 kDa dextran (which has a similar diffusion coefficient to VEGF) in the chamber is approximately constant after 15 h of diffusion in the central chamber of the experiment.

In this study, we develop a mathematical model based on the *in vitro* experimental results of cell migration in response to VEGF concentrations after 48 h. The endothelial cells were cultured for 24 h to enable endothelium formation before introducing VEGF into reservoirs on the other side of the central chamber, thereby creating an angiogenic factor gradient. Upon stimulation with an angiogenic factor gradient (0.1 *μg*/*mL*), a migration distance of 290 ± 53 *μm* was observed. This finding was derived in three separate microfluidic chips, each subjected to four measurements, resulting in a total of 12 migratory measurements. The experimental images can be found in Supplementary Figure S4 in the *“Microfluidic assay and cell culture preparation”* section of the Supplementary Material. In contrast, when no VEGF gradient was applied (0.0 *μg*/*mL*), the total cell migration was substantially lower, with a migration distance of 45 ± 26 *μm*. It is worth noting that the mean pore size of our hydrogel is 40 *μm*.The substantial increase in cell migration observed when the VEGF concentration was raised from 0.0 to 0.1 *μg*/*mL* highlights the profound impact of the VEGF gradient on endothelial cell migration and underscores the importance of VEGF in promoting angiogenesis. Based on these findings, our experiment indicates that stimulating cells with chemotactic factors is essential for effective migration. In the absence of a chemotactic factor gradient in the central chamber, cells do not migrate through the hydrogel but instead, fill the nearest pore.

In this study, the model focuses solely on a single cellular lineage, specifically endothelial cells, and does not account for other cell types, such as pericytes, which play a crucial role in vessel stabilisation and maturation. This simplification limits the applicability of the model to more complex *in vivo* situations. Second, the model only considers chemotaxis as the driving force behind cell migration without incorporating other forms of cell guidance, such as durotaxis or haptotaxis. The exclusion of these additional mechanisms may result in an incomplete understanding of the various factors contributing to cell migration *in vivo*.

The experiments conducted in this study aim to validate and compare results with those derived from our newly developed mathematical model. The chosen experimental conditions and mathematical model were specifically tailored to the system designed: we utilised a low concentration of fibrinogen, yielding a soft hydrogel with large pores conducive to cell size. It is important to note that this experimental setup represents only one aspect of the broader potential biological conditions, as it is intended to isolate the different migration contributors. In fact, we acknowledge that more complex, in vivo-like conditions, including varying sizes of pores, macropores, micropores, and different structures, might significantly impact cell migration patterns and angiogenesis.

## Model

This work proposes a dynamic phase-field model to describe the process of sprouting angiogenesis, where a collection of cells advances through a porous matrix responding to chemotactic factors. The model is based on a modified Cahn-Hilliard equation that incorporates the effects of the local environment on the angiogenic sprout growth. We are able to write in a single equation the whole migration model by only considering an extra contribution for the environmental effect without using agent-based terms. The governing equation is as follows:

Our formulation uses an order parameter, *ϕ*, to describe the presence or lack of cells in two phases taking values ±1 at the bulk. Here, *ϕ* represents the phase field variable characterising the angiogenic sprout interface, with *ϕ* = 1 denoting the sprout and *ϕ* = −1 indicating the surrounding environment. The term *ϕ* concerning time. *M* represents mobility, a proportionality constant that determines the speed of interface evolution. It can be thought of as a proportionality constant that relates the rate of change of the phase field variable to the driving forces of the system. The value of *ϵ* controls the width of the diffusive interface connecting the two phases, with a hyperbolic tangent profile, where *ϕ* values are close to zero. The first term of Eq. 1, the Cahn-Hilliard term, derives from a Ginzburg–Landau free energy written in dimensionless form Lacasta et al. (1993). The Cahn-Hilliard term, −*ϕ* + *ϕ*^{3} − *ϵ*^{2}∇^{2}*ϕ*, is responsible for the energy minimisation of the system. The terms −*ϕ* and *ϕ*^{3} respectively favor a homogeneous distribution and phase separation between the sprout and its environment. The gradient energy term, *ϵ*^{2}∇^{2}*ϕ*, penalizes large gradients of the phase field variable, promoting smooth interfaces. The environment effect term, 2*ϕϵC*_{0}*B*_{ϕ} − *ϵ*^{2}∇^{2}*B*_{ϕ}, accounts for the influence of the local environment on angiogenic sprout growth. *B*_{ϕ} denotes the biochemical signal guiding the sprouting process. The term 2*ϕϵC*_{0}*B*_{ϕ} represents the coupling between the phase field variable and the biochemical signal, promoting or inhibiting sprout growth depending on the local concentration of *B*_{ϕ}. The term −*ϵ*^{2}∇^{2}*B*_{ϕ} indicates the effect of the gradient of the biochemical signal on sprout growth, favouring growth in the direction of increasing *B*_{ϕ}. The total angiogenic factor, *C*_{0}, represents the environment in our model which incorporates the porous media and chemotactic factors, illustrated in Figure 2B), defined by

We assume that the dynamics of cells as a whole collection result from the balance between the biological agents present in their environment: the distribution of chemotactic factors (*ν*_{C}) and the presence of the extracellular matrix (*ν*_{ECM}). This reduced parameterisation improves the interpretability of the model by grouping the effects of the chemoattractants and the ECM. It is important to note that while the ECM can support cell migration by providing a physical framework, it can also hinder migration by creating barriers and obstacles for cells to overcome. Including the chemoattractants and the ECM in the model can capture this complexity and provide a more accurate representation of cell migration dynamics.

*C*_{0} represents three different environmental situations in our system: hydrogel alone, hydrogel combined with chemotactic factors and pore with chemotactic factors. In this article, if the contrary is not said, *ν*_{c} is taken as a linear gradient distribution (*ν*_{c}(*y*) = *αy* + *β*) and *ν*_{ECM} is based on the microscopy image of the hydrogel presented in (Figure 2A). Cells advance naturally if *C*_{0} is positive, but if it is zero or close to zero, cells do not feel its effect and remain immobile, allowing us to define the different situations cells encounter in the environment, (Figure 2B), quantitatively. The dark purple zones represent the densest ECM, or hydrogel, regions *ν*_{c} that activates cell migration

Following the description of Eq. 1, *B*_{ϕ} is the interaction term describing the influence of the cellular environment on cells. Tip cells lead the migration to proangiogenic and chemotactic cues present in their microenvironment Siekmann et al. (2013); Isogai et al. (2003); Gerhardt et al. (2003). Therefore, accounting for the minimal proliferation rate of tip cells Siemerink et al. (2012); Yetkin-Arik et al. (2019) and the time scale of our model, we have dismissed proliferation and allowed cells to migrate from the endothelium, modelled as a reservoir by using Dirichlet boundary conditions. *B*_{ϕ}, which communicates the effect of the microenvironment on the model, reads as:

The first term of the environment effect, 2*ϕϵC*_{0}*B*_{ϕ}, is derived from a chemotactic potential that drives and nurtures cells. The second term of the environment effect, *ϵ*^{2}∇^{2}*B*_{ϕ}, is a diffusive driving force that migrates cells at the cell interface. *B*_{ϕ} is proportional to *C*_{0} and centred at the limit of cells, keeping its effect far from the bulk cells, based on the rich cell-surface receptors and numerous filopodia present on tip cells. In our model, endothelial tip cells appear naturally by the localized impact of the biological factors at the interface; the non-homogeneity of the microenvironment promotes the advancement of cells with no additional rules. The values of the parameters used in the simulations are detailed in Table 1.

Our model conserves the number of cells in the system by using a conserved phase-field. Therefore, the initial condition determines the number of cells in the system and the number of mobile cells. However, to enable the advancement of cells as a collection, we consider the endothelium at the side-channel of the microfluidic platform acting as a reservoir of cells that injects cells into the system if necessary; cells stimulated chemotactically for advancement, *C*_{0} > 0 naturally come out from the reservoir. Dirichlet boundary conditions model this reservoir. To compute this migration flux, we use the total volume of the system at a given time in terms of its scalar field *ϕ* and divide it by the time:

Because the migration of ECs is oriented and regulated mechanically and chemically, we assume that the tip cell directs migration towards a clearer path. Moreover, EC collective migration also involves the degradation of the ECM to enable the progression of cells Lamalice et al. (2007); Yana et al. (2007); Koziol et al. (2012); Minor and Coulombe (2020); Song et al. (2018). For this reason, we extended the model by introducing a rule for tip ECs, allowing them to also degrade the ECM by producing matrix metalloproteinases (MMPs) and generating proangiogenic factors Seiki (2002); Haas and Madri (1999).

Several limitations are intrinsic to the mathematical model. Currently, the mathematical model is implemented in a two-dimensional (2D) framework, which, while simplifying the computational problem, has inherent limitations in its ability to represent the full complexity of 3D EC behaviours. We further assumed that the hydrogel is a homogeneous and isotropic material which does not accurately reflect the actual physical and mechanical properties of the fibrin-based hydrogel. Furthermore, the model considers only one type of cell. The interactions between endothelial cells and the hydrogel are simplified as chemoattractant gradients, neglecting mechanical interactions or adhesion between cells and the hydrogel. Lastly, the model is limited to a two-dimensional system and treats endothelial cells as a continuous phase without explicitly modelling individual cells, which may not accurately represent the three-dimensional complexity of the cellular microenvironment *in vivo*.

### Degradation

The ECM is a physical object directing, promoting and hindering cell motility depending on its physical and mechanical properties (porosity, pore size, rigidity, etc.) Zhang Y. et al. (2021); Cao et al. (2021); Bružauskaitė et al. (2016); Narayan and Venkatraman (2008); Qazi et al. (2019). Studies indicated that endothelial cells could sense patterns of extracellular distributions and interpret gradient directionality in 3D collagen gel for distances of the order of 100 *μm* Korff and Augustin (1999); Gerhardt et al. (2003). Moreover, the fact that the production of MMPs influences EC migration has been known for a long time. Therefore, based on the experimental observations, we also include in this model the ability of tip ECs to sense their environment and degrade the extracellular matrix through MMPs Hosseini et al. (2015); Gálvez et al. (2001). An example of how this mechanism is implemented is shown in Figure 3. Because of the bi-directional cell-ECM interaction, during the simulations of the extended model with degradation, the ECM is remodelled and re-biosynthesised as the tip ECs migrate Davis and Senger (2005). We assume that tip cells can sense the microenvironment approximately at a two-cell distance, approximately 40 *μm* from its tip, to direct migration. We base this behaviour on the combination of the multiple cell mechanisms combined to help acquire directionality and guide the cell pathfinding in response to the ECM cues (focal adhesions, filopodia, lamellipodia, dactylopodia, tip endothelial cells surface receptors and stalk cells receptors) Fischer et al. (2019); Figueiredo et al. (2021); Phng et al. (2013); Smet et al. (2009). Our model allow activated tip cells to avoid highly dense ECM regions and migrate through the ECM, taking advantage of the pores. More information about implementing this mechanism into the model can be found in the Supplementary Material.

**FIGURE 3**. Simulation of hydrogel testing and degradation for a static tip cell. The simulation shows snapshots of the ECM (purple) and chemotactic factor (green) at the tip sensing area (red square) for an immobile tip cell. The tip cell can sense approximately two cell distance (10-pixel distance) to encounter the path of lower ECM following the increasing chemotactic factor. It liberates matrix-degradation enzymes such as matrix metalloproteinases (MMPs), modelled in a Gaussian-like distribution, to increase the pore size. The first row shows the whole ECM domain, and the bottom row shows the tip cell sensing environment and its surrounding evolution when, without advancing, the tip cell produces MMPs. We show 0, 5 and 10 time steps from left to right for a 0.2 degradative power.

### Numerical integration

The implementation of our model is based on a 2D framework. This decision was made for simplification and computational efficiency, aiming to shed light on the complex dynamics of sprouting ECs. However, we emphasise that our 2D-implemented model inherently possesses limitations. Although it has provided valuable insights, the translation of these results to a more physiologically realistic three-dimensional (3D) environment should be approached with caution. In the 2D implemented model, cell behaviours and cell-cell and cell-matrix interactions might be simplified compared to those in a 3D setting. For example, cell migration patterns, as well as the distribution and perception of chemotactic gradients, could differ significantly between 2D and 3D contexts. Our ongoing research involves extending this model into the third dimension to capture more realistic cellular behaviours and interactions in a 3D environment.

In all simulations, reservoirs of the two phases have been placed by imposing Dirichlet boundary conditions. The reservoir of cells (*ϕ* = +1) has been placed at the boundary of the system for the simulations are shown in Figures 4, 5 and with the same distribution as the initial condition for simulations shown in Figures 6, 7.

**FIGURE 4**. Simplified pore system. The simulation in the top row, **(A)**, shows snapshots of the cells (red) far from the pore and chemotactic factor concentration (green to white). The middle row, **(B)**, shows the results of an analogous simulation in which the pore and biological factors are in contact with cells. Cells, *ϕ* = 1, advance from right to left, i.e., from *y* = 200 towards *y* =0. For **(A,B)**, we have the initial condition on the left, and on the right, we have the system at time *t* = 27.5*h*. Cells advance towards the maximum chemotactic factor concentration, i.e., from *y* = 200 towards *y* = 0. Red represents the position of cells (*ϕ* = +1) and blue (*ϕ* = −1) the non-cellular region. The bottom row, **(C)**, illustrates the position evolution of the tip of the cells in **(B)** at the centre of the vertical axis, i.e., *x* = 100. Cells respond only to local biological factors, advancing towards the maximum chemotactic factor concentration. The advancement of cells is fastest at initial times with a clear direction of the gradient of chemotactic concentration. In later stages, the velocity of advancement reduces, and the tip expands, occupying the pore region. Displayed in green is the chemotactic concentration distribution, and in blue is the tip position, both sharing the same *x*-axis *y*. We displayed the simulation time for the tip cell position on the left vertical axis. In the right vertical axis, we displayed the chemotactic factor concentration. The yellow arrow indicates the direction of advancement.

**FIGURE 5**. Fundamental dynamic model with lateral reservoir. **(A)** Represents the domain of the problem based on Figure 2) with an open path. *C*_{0} includes the effect of the extracellular matrix and the linear gradient of chemoattractant factors, increasing from right to left. **(B)** Cell advancement through the *C*_{0}-landscape shown in panel **(A)**, from the *C*_{0} minimum on the right to the *C*_{0} maximum on the left. The red area represents the initial condition of a spherical endothelium, and the red perimeter represents the enclosing cellular region at simulation time *t* = 96.25 h.

**FIGURE 6**. Fundamental dynamic model. Collective cell migration, with a small spherical reservoir and initial condition (filled in red), through an open path with no extracellular matrix degradation at *in silico* times {2.75, 27.75, 55, 96.25} hours, respectively from left to right. In red is the perimeter of cells advancing through the system.

**FIGURE 7**. Dynamic model with degradation. Collective cell migration with extracellular matrix degradation at time steps {2.75, 27.75, 55, 82.5} hours, respectively from left to right. We see the perimeter of cells in red advancing through the system. The collective cell migration can sense the densest regions of the extracellular matrix and avoid them while breaking down and migrating through lower-density extracellular matrix regions in the system.

To scale our *in silico* system to the represented physical system, we used the actual measures of the fluorescence microscopy image of the fibrin-based hydrogel. The hydrogel image has a height and width of 775*μm*, and we converted it into a 200*x*200 pixel image, giving a pixel length of 3.875 *μm*, which we truncated to 3.85 *μm*. Therefore our system size is *L* = 200 = 770 *μm*. The unit of time *in silico* is considered to be *t* = 3.0 ⋅ 10^{−3}*h*. All the rest of the parameters have been converted accordingly. *ϵ*, a parameter controlling the interface of our phase-field system, is set to be small compared to the size of the system (*ϵ* = 1 ≪ *L*) for proper convergence Rogers et al. (1988). If not stated otherwise, the Table 1 show all the parameter values used in our simulations. The continuum model is implemented in a vectorial form using Python and solved using the five-point stencil for spatial finite differences and an Euler scheme for the time dependence.We summarise in Table 1 the dimensionless values of the parameters used in the model and the corresponding physical values obtained using the size of system *L* = 770 *μm* and the total duration of the experiment, *T* = 48*h*. The mobility of the phase field that defines the location and advancement of endothelial cells is set by the timescale of the dynamics of capillaries, defined through Δ*t*.

## Results

(Figure 1A) depicts our experimental arrangement. As demonstrated in (Figure 1C), within our experimental system, the protrusive structures of migrating cells stretch to 290 ± 53 *μm* upon stimulation with a VEGF gradient of 0.1 *μg*/*mL* over 2 days. In contrast, when the system was not stimulated with an angiogenic factor gradient (0.0 *μg*/*mL*), cell migration was significantly reduced to a distance of 45 ± 26 *μm*. This sustained migration pattern gives rise to the capillary structures observed during angiogenesis. This pattern closely emulates natural capillaries and effectively demonstrates the characteristics of early-stage angiogenesis. In our in-silico model, we incorporated these experimental findings and developed a mathematical model to simulate cell migration in response to VEGF gradients. The model successfully replicated the experimentally observed migration distances under both conditions. In our migratory assays, there is some variation in cell migration through the extracellular matrix (ECM), as shown in (Figure 2A). This variation in collective cell migration can be attributed to the random structural differences present across various regions of the ECM, a factor that significantly influences cell migration.

To work with our model, we begin with a digitalised description of the cell migration results shown in (Figure 1C) by gray scaling and normalising the image and selecting a colourmap to highlight the presence of the cytoskeleton. Figures 1B, C) show cytoskeleton distribution for cells after 48 h without VEGF and under the influence of VEGF concentration gradient in our microfluidic assay. (Figure 1C) shows random small features indicating partial penetration and exploration of the cells into the ECM.

To gain mechanistic insights into various global and local behaviour of EC migration in our system, we have performed phase-field modelling. Here, we concentrate on the coupled influence of the ECM and the effect of VEGF, chemotactic agents, in cell migration. The ECM, by definition, is the surrounding material providing structural and biochemical support. The environment of cells in our phase-field model is comprised of two components: the hydrogel structure and the angiogenic factor distribution, *ν*_{ECM} and *ν*_{C}, respectively. In the following, we explore the advancement of cells reproduced by coupling these components into the *C*_{0} coefficient. We discover the apparent modulation effects of the ECM structure on the final cell morphology. Further, we allow cells to sense and degrade the ECM structure to reproduce cell migration.

Figure 4 presents a simplified system for cell advancement in a pore filled with chemotactic factors in the modelled ECM. In the first row, sprouting angiogenesis is not triggered due to the spatial distance between cells and the biological factors (the pore and the chemotactic factor). By contrast, panel B) in the second row shows a system where the pore filled up with chemotactic factor is in contact with cells, triggering invasion of cells into the pore region. Cells advance only with the local presence of positive values of total factor concentration, *C*_{0} (chemotactic factor and no dense ECM). The overall migration depends on the total angiogenic factor distribution (*C*_{0}) and the initial morphology of cells. Our model captures all essential attributes of cellular migration observed in our microfluidic assay, including pore and ECM invasion and chemotactic response. (Figure 4C) depicts the tip cell position (most advanced position *y*-direction wise) across the system and the *C*_{0} distribution as the chemotactic factor increases from 0.1, on the initial position of cells, to 0.5 and then decays to 0 again. Cells gradually invade the pore to follow the gradient direction of the angiogenic factor distribution. After the chemotactic distribution peak, the gradient direction disappears, resulting in a more isotropic cell advancement. The cell velocity varies with the change of *C*_{0} value and distribution. The velocity of cells is high until they reach the centre of the distribution, where their migration morphology changes from a sprouting penetrating tip towards an area-widening invasion to fill up the pore area. In this second regime, the advancing tip velocity substantially decreases because migrating cells no longer move towards the *y*-direction but homogeneously to the empty pore. A further decrease in cell migration velocity happens when the *C*_{0} value decreases below 0.1. The drop in cell velocity below a specific value indicates the lower bound of the total angiogenic factor, representing the interaction between cells and the ECM. Besides, the shape modulation of the ECM also causes cellular morphology adaptation, presumably through mechanical interactions of the actual system.

Figure 5 shows a numerical example of a representation of a domain of 750 *μm* by 750 *μm* of the central chamber in the experimental system. It has the same configuration as shown in (Figure 2B): a hydrogel structure with a general chemotactic factor gradient. In this case, the environment is static; cells are not able to modify it. For this reason, we have constructed a biomimetic pore in the ECM which cells follow. If this path was not created, the endothelium under the influence of this same environment as in Figure 2B), advances as a whole layer, activated by *ν*_{c}, following the chemotactic factor gradient. Activated migratory cells, unable to change their environment, sense the effect of the hydrogel generating small protrusions according to the ECM pores and irregularities. However, these protrusions are quickly absorbed by the advancing front of cells. However, by placing an open path, cells generate an invading protrusion that advances faster than those outside the pore that encounter the irregular ECM. Figure 5 shows that migrating cells adapt their shape to the morphology determined by the hydrogel while advancing towards the increasing chemotactic gradient. This result can be directly compared to Figure 1C, where cells migrated collectively towards a concentration gradient. Looking at the top and bottom, *x* = 0, 200, of Figure 5, we can also see that cells that do not encounter a clear open path also advance due to the chemotactic factor activation but at a much slower pace because of the hindering of the ECM.

This system, Figure 5, also explicitly shows how cells can advance over the irregular hydrogel pores slower than through the open trail. This tip cell advancement is comparable not only to Figure 1D but also to the results obtained by Wang et al. (2021) where ECM-cell interactions regulate the sprout diameter. We illustrate the dynamics of our model in Figure 6, where we changed the initial condition to an initial sphere, which will also act as our cell reservoir, to solely focus on the tip advancement that is the region of interest in our model. We see that this system with a lower pool of cells as the initial condition advances at the same speed and gives the same morphology as Figure 5. This result is not surprising as the advancing mechanism’s driving force is effectively the pore’s chemotactic factor distribution. The unique difference is the volume change, given that in Figure 5 we have the lateral cues also growing.

### Dynamic model with degradation

The advancement of cells through the ECM depends on open pores and the remodelling of the surrounding. Cells placed in a porous material must open their path through the ECM to migrate by exerting forces onto the ECM and degrading it, producing MMPs. For this reason, we extend our model by allowing tip cells to sense their environment and degrade, or re-biosynthesise, the hydrogel, converting it into pro-migratory agents Naito (2000). For simplicity, we have assumed that the production of MMPs from the tip is enough to avoid advancement hindered due to the presence of ECM. Using the capability of tip cells to sense their surroundings, we place the MMPs distribution according to the hydrogel concentration region sensed by the tip. Therefore the tip cell degrades the ECM, avoiding highly concentrated hydrogel regions as seen in Figure 7.

We can represent multiple scenarios by tailoring the MMPs production and its centring variability concerning the tip position. It is important to note that the degradative power and MMPs production play an essential role in deciding the tip path.

In Table 2 we summarize our experimental and numerical results. Consistently with our experiments, our model for cell migration predicts tip velocities of the same order of those obtained in our experiments. Because the driving mechanism for migration, chemotaxis, has the same profile in both models, we have obtained very similar mean velocities (differences due to path tortuosity) and migration rates (phase volume difference over time ratio) for the fundamental model and the model with degradation. In this model, the migration speed can be tailored to specific cases by changing the background gradient and the initial conditions of the system.

**TABLE 2**. System results. The *fundamental dynamic model* refers to Figure 5.A2-B2) and the *Dynamic model with degradation* refers to Figure 7. For each figure, average velocities have been computed between the initial and final conditions. The negative sign of mean tip y-velocity appears from cells advancing from *y* =200 towards *y* =0. *V* (*t*(*y*)) corresponds to the volume of cells in the system when the tip is at *x* as in Eq. 4 but subtracting the initial volume state and time, and *px* refers to pixel.

## Discussion

Our study does not consider the effects of fluid flow and shear stress on cell migration, as the design of our microfluidic setup aims to minimise these factors. However, it is important to recognise that under physiological conditions, both fluid flow and shear stress can significantly influence cell behaviour. Nonetheless, despite these constraints, our model offers crucial insights into the behaviour of endothelial cells in a controlled setting. These findings, in turn, can help inform future studies that may incorporate more nuanced and realistic microenvironment features.

To shed light on the fundamental interactions within an actual cell microenvironment, as well as to predict cell performance within a porous material, we developed a model grounded on microfluidic assays with rEPC, as detailed in Figure 1. The experimental model developed in this study applies to large-pore, low-concentration fibrin gels. It is here that the model’s robustness and applicability have been most effectively demonstrated. This platform enables us to isolate the effects of extracellular matrix (ECM) and the chemotactic factor, while effectively eliminating the influence of confounding parameters like fluid convection and flow Polacheck et al. (2011); Kingsmore et al. (2018).

The generation of capillaries is a complex process that highly depends on localised degradation of the ECM and subsequent EC migration, proliferation, and differentiation. ECs remodel the ECM during the invasion, and the ECM influences the shape and morphology of migratory cells. For ECs to form capillary-like structures, fibrin-based gels or ECM must be sufficiently soft to allow ECs to generate space to grow Kniazeva and Putnam (2009). The obtained results with our model elucidate the central role of the ECM during collective endothelial cell migration and organisation.

This work evaluated the capacity of the presented model to replicate the collective EC migration through a fibrin-based hydrogel. The hydrogel structure and composition, represented by *ν*_{ECM} in our model, are critical factors for the vascularisation process. In this work, we have focused on the structure of the hydrogel but have not compared or explored the adequacy of the specific ECM by changing the value of *ν*_{ECM} describing the ECM.

To introduce this new numerical model, we addressed the capability of chemotactic factors, such as VEGF, to promote cell migration in our experimental model. For this purpose, we tested our numerical system by testing the effect of a punctual diffusion source in a pore; see Figure 4. When ECs sense the chemotactic factor, they migrate towards the increasing gradient, agreeing with our experiments. On the other hand, for a homogeneous gradient, the cell advancement predicted in our model is also homogeneous, similar to epithelial sheet wound healing Nikolić et al. (2006). Introducing irregularities in this system, such as the ECM structure with no clear open paths or large pores, it is possible to further reproduce the small tips at the front of the epithelial sheet wound healing. If irregularities are not sustained during the advancement of cells, the advancing sheet of cells again absorbs the spontaneously formed tip structure. Therefore, to generate sprouting structures, we need a sustained inhomogeneity in our domain - which cells produce by degrading their surrounding. Similarly, we can use a punctual distribution of chemotactic factors or formed open paths through the ECM to describe sprouting-EC migration. Regions with irregular ECM hinder the advancement of ECs compared to open paths or pores, favouring fingering morphology led by the tip EC advancement as shown in Figures 5, 6. These finger-like shapes are a break of symmetry promoted by the surroundings. Therefore, bifurcations and multiple-fingering morphologies would naturally emerge given the appropriate surroundings with no need for extra rules conditioning the behaviour of the model.

The model with ECM degradation, shown in Figure 7, allows for dynamic matrix remodelling, an essential factor in angiogenesis. This process enables cells to form new pathways, leading to potentially more efficient and directed migration. The key distinction we observed relates to the tortuosity of the established pore, which impacts the *y*-axis tip velocity. By incorporating ECM degradation, we can see the effects of a variable and complex microenvironment on cell sprouting. However, the rate of degradation plays a crucial role. Overly rapid and extensive degradation could destroy the ECM, leading to unrealistic migration patterns. Conversely, a low degradation rate may restrict new pathway formation, thus limiting cell penetration. In contrast, the model without ECM degradation presumes a static ECM which takes advantage of the existing path. This limitation could lead to slower, less directed migration as cells are confined to the existing pore morphology.

The capability of cells to sense and migrate in their environment is a very complex process that influences cell migration. By allowing cells to sense and degrade their surrounding in our model with degradation, we can replicate the migratory behaviour seen in our experiments without manipulating the ECM structure, see Figure 7. Because we do not have more information in middle times, we set the degradation power of cells high enough not to hinder cell migration, as can be seen from the results shown in Table 2. At each time discretisation step, cells sense their surroundings, as in Figure 3, and decide to move accordingly to the path of lower environment resistance. We see how the average tip y-velocity in the case of ECM degradation is higher than the one with the pore already generated. The dynamic model with degradation opens a pore through the ECM more straight than the one defined for the dynamic model with no degradation. Besides, the phase field flux defined in Eq. 4 allows us to quantify the driving force stimulating the tip for a given set of initial conditions. The observed flux in both cases, with and without degradation, is of similar magnitude, albeit slightly higher in the case with ECM degradation. This outcome is anticipated since the chemotactic gradient profile remains the same for both scenarios. However, it is noteworthy that even a small portion of the ECM present in the path of migrating cells has the potential to contribute to subtle variations in cell sprouting behaviour and the establishment of gradients. This insight suggests that the degradation property and gradient profile can be readily adjusted to accommodate specific cases where the ECM may exhibit different tendencies toward degradation or remodelling.

While we have gleaned important insights into the dynamics of EC sprouting from our 2D implemented mathematical model, we recognise that these results have limitations when interpreted in a 3D physiological environment. For instance, cell behaviours and cell-cell and cell-matrix interactions in a 2D model can differ substantially from those in a more physiologically realistic 3D environment. A 3D environment might exhibit increased complexity in cell migration patterns due to the additional spatial degree of freedom, and the chemotactic gradients could distribute differently in a 3D context, influencing cell behaviour in ways not captured by our 2D implementation model.

We are actively working on extending our model into three dimensions to represent and predict EC behaviours in 3D experimental setups more accurately. Furthermore, we are employing additive manufacturing techniques together with different porogens to generate 3D structures with enhanced control and simulate a multitude of conditions with pores, in the micro and nanoscale. This will allow us in the future to compare the predictions of our mathematical model with cellular experiments under a wide range of conditions.

We have also confirmed that the dynamics or advancement of our model highly depends on the chemotactic factor distribution, which controls the flux of cells from the reservoir. Therefore, by scaling the chemotactic gradient and changing the cutoff of our ECM distribution, our model can reproduce multiple scenarios representing biomaterials with different morphologies. We believe that this model is an essential step toward modelling ECM-guided cell migration. Future developments and experiments should examine the effects of changes in the ECM on migrating cells and mechanical deformations on the same ECM structure.

## Conclusion

We proposed a phase-field model for cell migration that considers the ECM structure and chemotactic factors. From the mechanical point of view, the ECM morphology greatly influences the behaviour of cells by blocking and supporting their advancement. We compared our numerical model with experiments where endothelial cells migrate through the ECM from an endothelium formed in a side channel of the microfluidic system. In our model, we incorporated the ECM structure together with the chemotactic effect into the same parameter, describing the total angiogenic factor. With only one dynamic equation for our phase field, we can reproduce cellular migration and predict dynamics for the conformation of tip-like morphologies in the extracellular matrix. We also show that in our model, migrating cells adapt their shape to the morphology determined by the hydrogel while advancing towards the increasing chemotactic gradient. We have not included a proliferation term in our model; proliferation plays a significant role in collective cell dynamics at the lumen and capillary formation but not in endothelial tip advancement Siemerink et al. (2012); Yetkin-Arik et al. (2019).

We presented a fundamental dynamic model with a static environment and a dynamic model with extracellular matrix degradation for chemotactic cell migration. This novel model describes cell migration with a single equation and considers the extracellular matrix and chemotactic factors into a single parameter *C*_{0}. This new consideration reduces parametrization while increasing the interpretability of the model. In the first case, the fundamental dynamic model, where cells are not allowed to degrade the ECM, endothelial cells migrate through the pore, avoiding the ECM and producing a tip and fingering morphology. Multiple finger-like morphologies and bifurcations will naturally emerge if the environment presents suitable conditions without introducing any additional rules. By introducing degradation and micro-environment sensing in our dynamic model with degradation, our cells migrated towards the general chemotactic gradient avoiding highly dense ECM regions and degrading the necessary ECM to self-generate a pore. In all cases, the model predicted realistic dynamics.

Our model of endothelial cell migration predicts tip velocities of the same order as in the experiments. Additionally, it can be used to predict the behaviour of directed migrating cells through different extracellular matrices and chemotactic factor distributions by simply modifying the value of the coupling parameter, *C*_{0}. By including both the chemoattractants and the ECM in the model through this term, we can capture this complexity and provide a more accurate understanding of the dynamics of cell migration. This makes our model an excellent choice for testing and interpreting the behaviour of directed migrating cells in various microenvironments in laboratory experiments.

The model can be further improved, as this paper only shows the potential of the model to study the influence of the cellular environment on chemotactic cell migration. Our model opens a new framework to study cell migration interaction with different biomaterials. A simple way to extend the model would be to consider deformable extracellular matrices through a viscoelastic model to account for the rheology and deformations of the structure, leading to a more realistic dynamics representation.

## Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

## Ethics statement

All Bone marrow-derived rat endothelial progenitor cells (rEPCs) used in this study were properly obtained from long bones of young Lewis rats (2-4 weeks old) following methods in accordance with relevant guidelines and regulations, and all animal care protocols were approved by the Committee on Ethics and Animal Experiments of the Scientific Park of Barcelona (Permit No. 0006S/13393/2011, 2011) Aguirre et al. (2010); Gonzalez-Vazquez et al. (2014).

## Author contributions

JF-T performed all the numerical simulations. AN-M performed the experiments. AL-C designed the microfluidic platform. OC designed the experiments. JF-T, JR-A, RB, and AH-M developed the model. JR-A, RB, OC, and AH-M revised the project’s main conceptual ideas. All authors contributed to the article and approved the submitted version.

## Funding

This work was supported in part by the Spanish Ministry of Science and Innovation (MICINN) and the State Research Agency (AEI) through the projects RTI 2018-097038-B-C22, PID 2019-106063GB-I00, PID2021-124575OB-I00 and CONACyT through the project 283279. JF-T is supported by Generalitat de Catalunya/AGAUR 2018-DI-064. AL-C thanks The Spanish Ministry of Universities for the fellowship FPU17/06161. JR-A acknowledges support from DGAPA-UNAM grant IA-100823. AH-M is supported by Generalitat de Catalunya/AGAUR 2021SGR00450. Authors also thank for the grant from the University of Barcelona to publish in open access.

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2023.1145550/full#supplementary-material

## References

Aguirre, A., González, A., Planell, J., and Engel, E. (2010). Extracellular calcium modulates *in vitro* bone marrow-derived flk-1+ cd34+ progenitor cell chemotaxis and differentiation through a calcium-sensing receptor. *Biochem. biophysical Res. Commun.* 393, 156–161. doi:10.1016/j.bbrc.2010.01.109

Alert, R., and Trepat, X. (2020). Physical models of collective cell migration. *Annu. Rev. Condens. Matter Phys.* 11, 77–101. doi:10.1146/annurev-conmatphys-031218-013516

Blache, U., and Ehrbar, M. (2018). Inspired by nature: Hydrogels as versatile tools for vascular engineering. *Adv. wound care* 7, 232–246. doi:10.1089/wound.2017.0760

Blocki, A., Beyer, S., Jung, F., and Raghunath, M. (2018). The controversial origin of pericytes during angiogenesis–implications for cell-based therapeutic angiogenesis and cell-based therapies. *Clin. Hemorheol. Microcirc.* 69, 215–232. doi:10.3233/ch-189132

Bružauskaitė, I., Bironaitė, D., Bagdonas, E., and Bernotienė, E. (2016). Scaffolds and cells for tissue regeneration: Different scaffold pore sizes—different cell effects. *Cytotechnology* 68, 355–369. doi:10.1007/s10616-015-9895-4

Caló, E., and Khutoryanskiy, V. V. (2015). Biomedical applications of hydrogels: A review of patents and commercial products. *Eur. Polym. J.* 65, 252–267. doi:10.1016/j.eurpolymj.2014.11.024

Cao, Y., Jiang, J., Jiang, Y., Li, Z., Hou, J., and Li, Q. (2021). Biodegradable highly porous interconnected poly(*ɛ*-caprolactone)/poly(l-lactide-co-*ɛ*-caprolactone) scaffolds by supercritical foaming for small-diameter vascular tissue engineering. *Polym. Adv. Technol.* 33, 440–451. doi:10.1002/pat.5528

Carmeliet, P., De Smet, F., Loges, S., and Mazzone, M. (2009). Branching morphogenesis and antiangiogenesis candidates: Tip cells lead the way. *Nat. Rev. Clin. Oncol.* 6, 315–326. doi:10.1038/nrclinonc.2009.64

Carmeliet, P., and Jain, R. K. (2011). Molecular mechanisms and clinical applications of angiogenesis. *Nature* 473, 298–307. doi:10.1038/nature10144

Davis, G. E., and Senger, D. R. (2005). Endothelial extracellular matrix. *Circulation Res.* 97, 1093–1107. doi:10.1161/01.RES.0000191547.64391.e3

De Smet, F., Segura, I., De Bock, K., Hohensinner, P. J., and Carmeliet, P. (2009). Mechanisms of vessel branching: Filopodia on endothelial tip cells lead the way. *Arteriosclerosis, thrombosis, Vasc. Biol.* 29, 639–649. doi:10.1161/atvbaha.109.185165

Elosegui-Artola, A. (2021). The extracellular matrix viscoelasticity as a regulator of cell and tissue dynamics. *Curr. Opin. Cell Biol.* 72, 10–18. doi:10.1016/j.ceb.2021.04.002

Figueiredo, A. M., Barbacena, P., Russo, A., Vaccaro, S., Ramalho, D., Pena, A., et al. (2021). Endothelial cell invasion is controlled by dactylopodia. *Proc. Natl. Acad. Sci.* 118, e2023829118. doi:10.1073/pnas.2023829118

Fischer, R. S., Lam, P.-Y., Huttenlocher, A., and Waterman, C. M. (2019). Filopodia and focal adhesions: An integrated system driving branching morphogenesis in neuronal pathfinding and angiogenesis. *Dev. Biol.* 451, 86–95. doi:10.1016/j.ydbio.2018.08.015

Friedl, P., and Gilmour, D. (2009). Collective cell migration in morphogenesis, regeneration and cancer. *Nat. Rev. Mol. cell Biol.* 10, 445–457. doi:10.1038/nrm2720

Gálvez, B. G., Matiás-Román, S., Albar, J. P., Sánchez-Madrid, F., and Arroyo, A. G. (2001). Membrane type 1-matrix metalloproteinase is activated during migration of human endothelial cells and modulates endothelial motility and matrix remodeling. *J. Biol. Chem.* 276, 37491–37500. doi:10.1074/jbc.m104094200

Gerhardt, H., Golding, M., Fruttiger, M., Ruhrberg, C., Lundkvist, A., Abramsson, A., et al. (2003). VEGF guides angiogenic sprouting utilizing endothelial tip cell filopodia. *J. Cell Biol.* 161, 1163–1177. doi:10.1083/jcb.200302047

Ghaffarizadeh, A., Heiland, R., Friedman, S. H., Mumenthaler, S. M., and Macklin, P. (2018). Physicell: An open source physics-based cell simulator for 3-d multicellular systems. *PLoS Comput. Biol.* 14, e1005991. doi:10.1371/journal.pcbi.1005991

Giverso, C., and Ciarletta, P. (2016). Erratum: Tumour angiogenesis as a chemo-mechanical surface instability. *Sci. Rep.* 6, 25214. doi:10.1038/srep25214

González-Vázquez, A., Planell, J. A., and Engel, E. (2014). Extracellular calcium and casr drive osteoinduction in mesenchymal stromal cells. *Acta biomater.* 10, 2824–2833. doi:10.1016/j.actbio.2014.02.004

Graner, F., and Glazier, J. A. (1992). Simulation of biological cell sorting using a two-dimensional extended potts model. *Phys. Rev. Lett.* 69, 2013–2016. doi:10.1103/physrevlett.69.2013

Haas, T., and Madri, J. A. (1999). Extracellular matrix-driven matrix metalloproteinase production in endothelial cells: Implications for angiogenesis. *Trends Cardiovasc. Med.* 9, 70–77. doi:10.1016/s1050-1738(99)00014-6

Heck, T., Vaeyens, M.-M., and Van Oosterwyck, H. (2015). Computational models of sprouting angiogenesis and cell migration: Towards multiscale mechanochemical models of angiogenesis. *Math. Model. Nat. Phenom.* 10, 108–141. doi:10.1051/mmnp/201510106

Hendrata, M., and Sudiono, J. (2019). A hybrid multiscale model for investigating tumor angiogenesis and its response to cell-based therapy. *Silico Biol.* 13, 1–20. doi:10.3233/isb-170469

Hosseini, Y., Agah, M., and Verbridge, S. S. (2015). Endothelial cell sensing, restructuring, and invasion in collagen hydrogel structures. *Integr. Biol.* 7, 1432–1441. doi:10.1039/c5ib00207a

Isogai, S., Lawson, N. D., Torrealday, S., Horiguchi, M., and Weinstein, B. M. (2003). Angiogenic network formation in the developing vertebrate trunk. *Development* 130, 5281–5290. doi:10.1242/dev.00733

Jain, R. K. (2003). Molecular regulation of vessel maturation. *Nat. Med.* 9, 685–693. doi:10.1038/nm0603-685

Jakobsson, L., Franco, C. A., Bentley, K., Collins, R. T., Ponsioen, B., Aspalter, I. M., et al. (2010). Endothelial cells dynamically compete for the tip cell position during angiogenic sprouting. *Nat. cell Biol.* 12, 943–953. doi:10.1038/ncb2103

Kanchanawong, P., and Calderwood, D. A. (2023). Organization, dynamics and mechanoregulation of integrin-mediated cell–ecm adhesions. *Nat. Rev. Mol. Cell Biol.* 24, 142–161. doi:10.1038/s41580-022-00531-5

Kingsmore, K. M., Vaccari, A., Abler, D., Cui, S. X., Epstein, F. H., Rockne, R. C., et al. (2018). Mri analysis to map interstitial flow in the brain tumor microenvironment. *Apl. Bioeng.* 2, 031905. doi:10.1063/1.5023503

Kniazeva, E., and Putnam, A. J. (2009). Endothelial cell traction and ecm density influence both capillary morphogenesis and maintenance in 3-d. *Am. J. Physiology-Cell Physiology* 297, C179–C187. doi:10.1152/ajpcell.00018.2009

Korff, T., and Augustin, H. G. (1999). Tensional forces in fibrillar extracellular matrices control directional capillary sprouting. *J. cell Sci.* 112, 3249–3258. doi:10.1242/jcs.112.19.3249

Koziol, A., Gonzalo, P., Mota, A., Pollán, Á., Lorenzo, C., Colomé, N., et al. (2012). The protease mt1-mmp drives a combinatorial proteolytic program in activated endothelial cells. *FASEB J.* 26, 4481–4494. doi:10.1096/fj.12-205906

Krishnamoorthy, S., Zhang, Z., and Xu, C. (2019). Biofabrication of three-dimensional cellular structures based on gelatin methacrylate–alginate interpenetrating network hydrogel. *J. Biomaterials Appl.* 33, 1105–1117. doi:10.1177/0885328218823329

Krüger-Genge, A., Blocki, A., Franke, R.-P., and Jung, F. (2019). Vascular endothelial cell biology: An update. *Int. J. Mol. Sci.* 20, 4411. doi:10.3390/ijms20184411

Lacasta, A., Hernández-Machado, A., and Sancho, J. M. (1993). Front and domain growth in the presence of gravity. *Phys. Rev. B* 48, 9418–9425. doi:10.1103/physrevb.48.9418

Lamalice, L., Boeuf, F. L., and Huot, J. (2007). Endothelial cell migration during angiogenesis. *Circulation Res.* 100, 782–794. doi:10.1161/01.RES.0000259593.07661.1e

Leschiera, E., Lorenzi, T., Shen, S., Almeida, L., and Audebert, C. (2022). A mathematical model to study the impact of intra-tumour heterogeneity on anti-tumour cd8+ t cell immune response. *J. Theor. Biol.* 538, 111028. doi:10.1016/j.jtbi.2022.111028

Lima, E., Oden, J., and Almeida, R. (2014). A hybrid ten-species phase-field model of tumor growth. *Math. Models Methods Appl. Sci.* 24, 2569–2599. doi:10.1142/s0218202514500304

López-Canosa, A., Pérez-Amodio, S., Engel, E., and Castaño, O. (2022). Microfluidic 3d platform to evaluate endothelial progenitor cell recruitment by bioactive materials. *Acta Biomater.* 151, 264–277. doi:10.1016/j.actbio.2022.08.019

López-Canosa, A., Perez-Amodio, S., Yanac-Huertas, E., Ordoño, J., Rodriguez-Trujillo, R., Samitier, J., et al. (2021). A microphysiological system combining electrospun fibers and electrical stimulation for the maturation of highly anisotropic cardiac tissue. *Biofabrication* 13, 035047. doi:10.1088/1758-5090/abff12

Lutolf, M., and Hubbell, J. (2005). Synthetic biomaterials as instructive extracellular microenvironments for morphogenesis in tissue engineering. *Nat. Biotechnol.* 23, 47–55. doi:10.1038/nbt1055

Ma, Y., Han, T., Yang, Q., Wang, J., Feng, B., Jia, Y., et al. (2021). Viscoelastic cell microenvironment: Hydrogel-based strategy for recapitulating dynamic ecm mechanics. *Adv. Funct. Mater.* 31, 2100848. doi:10.1002/adfm.202100848

Marchand, M., Monnot, C., Muller, L., and Germain, S. (2019). Extracellular matrix scaffolding in angiogenesis and capillary homeostasis. *Seminars cell and Dev. Biol. (Elsevier)* 89, 147–156. doi:10.1016/j.semcdb.2018.08.007

Mayor, R., and Etienne-Manneville, S. (2016). The front and rear of collective cell migration. *Nat. Rev. Mol. cell Biol.* 17, 97–109. doi:10.1038/nrm.2015.14

Minor, A. J., and Coulombe, K. L. (2020). Engineering a collagen matrix for cell-instructive regenerative angiogenesis. *J. Biomed. Mater. Res. Part B Appl. Biomaterials* 108, 2407–2416. doi:10.1002/jbm.b.34573

Moure, A., and Gomez, H. (2021). Phase-field modeling of individual and collective cell migration. *Archives Comput. Methods Eng.* 28, 311–344. doi:10.1007/s11831-019-09377-1

Naito, M. (2000). Effects of fibrinogen, fibrin and their degradation products on the behaviour of vascular smooth muscle cells. *Nihon Ronen Igakkai zasshi. Jpn. J. Geriatrics* 37, 458–463. doi:10.3143/geriatrics.37.458

Narayan, D., and Venkatraman, S. S. (2008). Effect of pore size and interpore distance on endothelial cell growth on polymers. *J. Biomed. Mater. Res. Part A* 87A, 710–718. doi:10.1002/jbm.a.31749

Nikolić, D. L., Boettiger, A. N., Bar-Sagi, D., Carbeck, J. D., and Shvartsman, S. Y. (2006). Role of boundary conditions in an experimental model of epithelial wound healing. *Am. J. Physiology-Cell Physiology* 291, C68–C75. doi:10.1152/ajpcell.00411.2005

Norton, K.-A., Gong, C., Jamalian, S., and Popel, A. S. (2019). Multiscale agent-based and hybrid modeling of the tumor immune microenvironment. *Processes* 7, 37. doi:10.3390/pr7010037

Phng, L.-K., Stanchi, F., and Gerhardt, H. (2013). Filopodia are dispensable for endothelial tip cell guidance. *Development* 140, 4031–4040. doi:10.1242/dev.097352

Pitulescu, M. E., Schmidt, I., Giaimo, B. D., Antoine, T., Berkenfeld, F., Ferrante, F., et al. (2017). Dll4 and notch signalling couples sprouting angiogenesis and artery formation. *Nat. cell Biol.* 19, 915–927. doi:10.1038/ncb3555

Polacheck, W. J., Charest, J. L., and Kamm, R. D. (2011). Interstitial flow influences direction of tumor cell migration through competing mechanisms. *Proc. Natl. Acad. Sci.* 108, 11115–11120. doi:10.1073/pnas.1103581108

Qazi, T. H., Tytgat, L., Dubruel, P., Duda, G. N., Van Vlierberghe, S., and Geissler, S. (2019). Extrusion printed scaffolds with varying pore size as modulators of msc angiogenic paracrine effects. *ACS Biomaterials Sci. Eng.* 5, 5348–5358. doi:10.1021/acsbiomaterials.9b00843

Qiu, J., and Hirschi, K. K. (2019). Endothelial cell development and its application to regenerative medicine. *Circulation Res.* 125, 489–501. doi:10.1161/circresaha.119.311405

Rogers, T., Elder, K., and Desai, R. C. (1988). Numerical study of the late stages of spinodal decomposition. *Phys. Rev. B* 37, 9638–9649. doi:10.1103/physrevb.37.9638

Santos-Oliveira, P., Correia, A., Rodrigues, T., Ribeiro-Rodrigues, T. M., Matafome, P., Rodríguez-Manzaneque, J. C., et al. (2015). The force at the tip-modelling tension and proliferation in sprouting angiogenesis. *PLoS Comput. Biol.* 11, e1004436. doi:10.1371/journal.pcbi.1004436

Santra, S., Mandal, S., and Chakraborty, S. (2020). Phase-field modeling of multicomponent and multiphase flows in microfluidic systems: A review. *Int. J. Numer. Methods Heat Fluid Flow* 31, 3089–3131. doi:10.1108/hff-01-2020-0001

Schultz, G. S., and Wysocki, A. (2009). Interactions between extracellular matrix and growth factors in wound healing. *Wound repair Regen.* 17, 153–162. doi:10.1111/j.1524-475x.2009.00466.x

Seiki, M. (2002). The cell surface: The stage for matrix metalloproteinase regulation of migration. *Curr. Opin. cell Biol.* 14, 624–632. doi:10.1016/s0955-0674(02)00363-0

Shibuya, M. (2001). Structure and function of vegf/vegf-receptor system involved in angiogenesis. *Cell Struct. Funct.* 26, 25–35. doi:10.1247/csf.26.25

Siekmann, A. F., Affolter, M., and Belting, H.-G. (2013). The tip cell concept 10 years after: New players tune in for a common theme. *Exp. cell Res.* 319, 1255–1263. doi:10.1016/j.yexcr.2013.01.019

Siemerink, M. J., Klaassen, I., Vogels, I. M., Griffioen, A. W., Van Noorden, C. J., and Schlingemann, R. O. (2012). Cd34 marks angiogenic tip cells in human vascular endothelial cell cultures. *Angiogenesis* 15, 151–163. doi:10.1007/s10456-011-9251-z

Smet, F. D., Segura, I., Bock, K. D., Hohensinner, P. J., and Carmeliet, P. (2009). Mechanisms of vessel branching. *Arteriosclerosis, Thrombosis, Vasc. Biol.* 29, 639–649. doi:10.1161/ATVBAHA.109.185165

Song, K. H., Highley, C. B., Rouff, A., and Burdick, J. A. (2018). Complex 3d-printed microchannels within cell-degradable hydrogels. *Adv. Funct. Mater.* 28, 1801331. doi:10.1002/adfm.201801331

Szabó, A., and Merks, R. M. (2013). Cellular potts modeling of tumor growth, tumor invasion, and tumor evolution. *Front. Oncol.* 3, 87. doi:10.3389/fonc.2013.00087

Travasso, R. D., Corvera Poiré, E., Castro, M., Rodrguez-Manzaneque, J. C., and Hernández-Machado, A. (2011). Tumor angiogenesis and vascular patterning: A mathematical model. *PloS one* 6, e19989. doi:10.1371/journal.pone.0019989

van Oers, R. F., Rens, E. G., LaValley, D. J., Reinhart-King, C. A., and Merks, R. M. (2014). Mechanical cell-matrix feedback explains pairwise and collective endothelial cell behavior *in vitro*. *PLoS Comput. Biol.* 10, e1003774. doi:10.1371/journal.pcbi.1003774

Vats, K., and Benoit, D. S. (2013). Dynamic manipulation of hydrogels to control cell behavior: A review. *Tissue Eng. Part B Rev.* 19, 455–469. doi:10.1089/ten.teb.2012.0716

Vila, A., Torras, N., Castaño, A. G., García-Díaz, M., Comelles, J., Pérez-Berezo, T., et al. (2020). Hydrogel co-networks of gelatine methacrylate and poly (ethylene glycol) diacrylate sustain 3d functional *in vitro* models of intestinal mucosa. *Biofabrication* 12, 025008. doi:10.1088/1758-5090/ab5f50

Vilanova, G., Colominas, I., and Gomez, H. (2017). A mathematical model of tumour angiogenesis: Growth, regression and regrowth. *J. R. Soc. Interface* 14, 20160918. doi:10.1098/rsif.2016.0918

Vilanova, G., Colominas, I., and Gomez, H. (2013). Capillary networks in tumor angiogenesis: From discrete endothelial cells to phase-field averaged descriptions via isogeometric analysis. *Int. J. Numer. methods Biomed. Eng.* 29, 1015–1037. doi:10.1002/cnm.2552

Wang, W. Y., Jarman, E. H., Lin, D., and Baker, B. M. (2021). Dynamic endothelial stalk cell–matrix interactions regulate angiogenic sprout diameter. *Front. Bioeng. Biotechnol.* 9, 620128. doi:10.3389/fbioe.2021.620128

Xu, J., Vilanova, G., and Gomez, H. (2016). A mathematical model coupling tumor growth and angiogenesis. *PloS one* 11, e0149422. doi:10.1371/journal.pone.0149422

Xu, J., Vilanova, G., and Gomez, H. (2017). Full-scale, three-dimensional simulation of early-stage tumor growth: The onset of malignancy. *Comput. Methods Appl. Mech. Eng.* 314, 126–146. doi:10.1016/j.cma.2016.07.010

Xu, J., Vilanova, G., and Gomez, H. (2020). Phase-field model of vascular tumor growth: Three-dimensional geometry of the vascular network and integration with imaging data. *Comput. Methods Appl. Mech. Eng.* 359, 112648. doi:10.1016/j.cma.2019.112648

Yana, I., Sagara, H., Takaki, S., Takatsu, K., Nakamura, K., Nakao, K., et al. (2007). Crosstalk between neovessels and mural cells directs the site-specific expression of MT1-MMP to endothelial tip cells. *J. Cell Sci.* 120, 1607–1614. doi:10.1242/jcs.000679

Yetkin-Arik, B., Vogels, I., Neyazi, N., Van Duinen, V., Houtkooper, R., Van Noorden, C., et al. (2019). Endothelial tip cells *in vitro* are less glycolytic and have a more flexible response to metabolic stress than non-tip cells. *Sci. Rep.* 9, 10414–10417. doi:10.1038/s41598-019-46503-2

Zhang, K., Feng, Q., Fang, Z., Gu, L., and Bian, L. (2021a). Structurally dynamic hydrogels for biomedical applications: Pursuing a fine balance between macroscopic stability and microscopic dynamics. *Chem. Rev.* 121, 11149–11193. doi:10.1021/acs.chemrev.1c00071

Zhang, Y., Wang, K., Du, Y., Yang, H., Jia, G., Huang, D., et al. (2023). Computational modeling to determine the effect of phenotypic heterogeneity in tumors on the collective tumor–immune interactions. *Bull. Math. Biol.* 85, 51. doi:10.1007/s11538-023-01158-z

Keywords: extracellular matrix, angiogenesis, chemotaxis, endothelial cells, biomimmetic, phase field, mathematical models, *in silico* model

Citation: Ferre-Torres J, Noguera-Monteagudo A, Lopez-Canosa A, Romero-Arias JR, Barrio R, Castaño O and Hernandez-Machado A (2023) Modelling of chemotactic sprouting endothelial cells through an extracellular matrix. *Front. Bioeng. Biotechnol.* 11:1145550. doi: 10.3389/fbioe.2023.1145550

Received: 17 January 2023; Accepted: 26 May 2023;

Published: 08 June 2023.

Edited by:

Arturo Ibáñez-Fonseca, Lund University, SwedenReviewed by:

Mathieu Hautefeuille, Sorbonne Université, FranceJoao Carlos Carvalho, University of Coimbra, Portugal

Copyright © 2023 Ferre-Torres, Noguera-Monteagudo, Lopez-Canosa, Romero-Arias, Barrio, Castaño and Hernandez-Machado. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Josep Ferre-Torres, josep.ferre@fmc.ub.edu; Aurora Hernandez-Machado, a.hernandezmachado@ub.edu