Multi-Cue Kinetic Model with Non-Local Sensing for Cell Migration on a Fiber Network with Chemotaxis.

Cells perform directed motion in response to external stimuli that they detect by sensing the environment with their membrane protrusions. Precisely, several biochemical and biophysical cues give rise to tactic migration in the direction of their specific targets. Thus, this defines a multi-cue environment in which cells have to sort and combine different, and potentially competitive, stimuli. We propose a non-local kinetic model for cell migration in which cell polarization is influenced simultaneously by two external factors: contact guidance and chemotaxis. We propose two different sensing strategies, and we analyze the two resulting transport kinetic models by recovering the appropriate macroscopic limit in different regimes, in order to observe how the cell size, with respect to the variation of both external fields, influences the overall behavior. This analysis shows the importance of dealing with hyperbolic models, rather than drift-diffusion ones. Moreover, we numerically integrate the kinetic transport equations in a two-dimensional setting in order to investigate qualitatively various scenarios. Finally, we show how our setting is able to reproduce some experimental results concerning the influence of topographical and chemical cues in directing cell motility.

Introduction

Cell migration is a fundamental mechanism in a huge variety of processes, such as wound healing, angiogenesis, tumor stroma formation, and metastasis. During these processes, cells sense the environment and respond to external stimuli orienting their direction of motion toward specific targets. This mechanism is referred to as taxis, and it results in the persistent migration in a certain preferential direction. The guidance cues leading to directed migration may be biochemical or biophysical. One example of a biochemical cue is the concentration of soluble molecules in the extracellular space. This cue gives rise to chemotaxis, which is considered a mono-directional stimulus. Other cues generating mono-directional stimuli include electric fields (electrotaxis, or galvanotaxis), light signals (phototaxis), bound ligands to the substratum (haptotaxis), or the extracellular matrix (ECM) stiffness (durotaxis) (Lara and Schneider 2013). Precisely, ECM stiffness can be counted as a biophysical cue, as well as the collagen fiber alignment. In particular, the latter is shown to stimulate contact guidance (Friedl and Brocker 2000; Friedl 2004), i.e., the tendency of cells to migrate by crawling on the fibers and following the directions imposed by them. Contact guidance is a bi-directional cue. In fact, if the fibers are not polarized, there is no preferential sense of motion along them.

In many pathological and physiological processes, there are several directional cues inducing different simultaneous stimuli. While the cell’s response to each of them has been largely studied, the cell’s response to a multi-cue environment is much less understood. Some studies have shown how there can be competition or cooperation between these different stimuli. Thus, the fundamental issue concerns the way cells rank, integrate, or hierarchize them, especially when these stimuli are competing (e.g., when they are not co-aligned) (Rajnicek et al. 2007). Therefore, with the present work we propose a kinetic model aimed at analyzing cell behavior in response to two different stimuli. We study the way the simultaneous sensing of two cues—chemotaxis and contact guidance—influences the choice of the cell migratory direction. We take into account non-local sensing of both cues, since cells extend their protrusions in order to sense the environmental stimuli.

Biological Background

The coexistence of chemotaxis and contact guidance happens in vivo in a variety of situations, like wound healing or cancer progression. For example, during wound healing, fibroblasts migrate efficiently along collagen or fibronectin fibers in connective tissues. In tumor spread and metastasis formation, cancer cells follow the aligned fibers at the tumor-stroma interface and, thus, are facilitated to reach blood and lymphatic vessels (Steeg 2016; Provenzano et al. 2006, 2009). In both cases, chemotactic gradients have been shown to accelerate and enhance these processes (Lara and Schneider 2013; Bromberek et al. 2002). Another important issue concerns the design of platforms for controlling multiple directional cues and, in particular, soluble factors and aligned fibers. In fact, there are not many experimental studies that look at the combined effect of chemotaxis and contact guidance (Lara and Schneider 2013). In one of the first works on this topic (Wilkinson and Lackie 1983), the authors analyze contact guidance of neutrophil leukocytes on fibrils of collagen, showing a more efficient migration in the fiber direction, instead of in the perpendicular one. They also observe that, in the presence of a chemoattractant, there is cooperation or competition between the cues depending on their relative orientations. In the work by Bromberek et al. (2002), the enhancement of the alignment along the fibers is observed in presence of a co-aligned chemoattractant, while, in Maheshwari et al. (1999), the authors study the effects of different fibronectin densities and growth factor (EGF) concentrations on the quantitative regulation of random cell migration. An interesting 2D platform, allowing to study contact guidance and chemotaxis, was proposed by Sundararaghavan et al. (2013). Here, the authors demonstrate an additive effect of chemical gradients and fiber alignment by measuring the persistence time and a stronger dominance of contact guidance when the chemotactic gradient is aligned perpendicularly to the fibers. Thus, as for multiple directional cues different scenarios may happen, a deep understanding of cell migrational responses is a key step for the comprehension of both physiological and pathological processes.

Mathematical Models for Directed Cell Migration

There is a huge variety of mathematical models concerning different aspects of directed cell migration. They range from microscopic (also called individual-based) models, describing migration at the cell level, up to macroscopic ones, concerning collective cell migration at a tissue level.

Among individual-based models, there are many examples regarding separately chemotaxis (e.g., see Di Costanzo et al. 2020; Giniūnaitė et al. 2019 and references therein) and migration on the ECM (Colombi et al. 2017; Scianna et al. 2013; Schlüter et al. 2012). Concerning macroscopic settings, the famous Keller–Segel model is one of the first examples of a drift-diffusion system postulated at the macroscopic level for the study of the chemotactic effects (Keller and Segel 1970). Then, many efforts were made in order to encompass some defects of this setting, as well as for deriving it from lower scale models (e.g., see Painter 2019; Hillen and Painter 2008; Othmer and Hillen 2002; Othmer and Stevens 2001; Bellomo et al. 2015 and references therein). Between microscopic and macroscopic models, the mesoscopic settings represent an intermediate scale, which includes microscopic dynamics and describes the statistical distribution of the individuals. As in the case of kinetic theory, these models allow to recover the appropriate macroscopic regime, which inherits some details of the microscopic dynamics and, thus, gives more significance to some of the derived parameters (e.g., see Othmer and Hillen 2002; Colombi et al. 2015; Chalub et al. 2004; Eftimie 2012 and references therein). The two major models describing contact guidance at the mesoscopic level were proposed by Dickinson (2000) and Hillen (2006), and both models are local in the physical space.

Although the wide literature concerning single-cue models, not many models exist for multiple cues, especially settings describing in details the process of cell sensing. Among the few examples of multiple-cue analysis, we recall the work proposed by Painter et al. (2000), where a macroscopic drift-diffusion model is derived from a space jump process describing the response to multiple chemicals. A review of macroscopic PDEs including multiple-taxis has been recently proposed by Kolbe et al. (2021). Looking specifically at the combination of chemotaxis and contact guidance, in Wagle and Tranquillo (2000) one of the first models describing microscopic dynamics is proposed, while, in Azimzade et al. (2019), a microscopic double-cue stochastic model was introduced to analyze cell migration and classify tumor associated collagen signatures (TACS). We also recall the kinetic models for cell–cell interactions on a fiber network in presence of a tactic cue presented by Chauviere et al. (2007). Recently, an extension of the M5 model (Hillen 2006) was extended in order to include chemotaxis, with the aim of showing network formation in the presence of a chemoattractant (Thiessen and Hillen 2021). Non-locality in the description of cell sensing over a finite neighborhood and in presence of two external cues was firstly introduced by Loy and Preziosi (2020, 2020). This non-local kinetic model takes into account the dependence of the cell sensing on the cell size (e.g., on the maximum cell protrusion length) and includes two cues, one affecting cell direction, while another influencing cell speed.

Outline of the Paper

Due to the lack of mathematical models looking at the interplay between chemotaxis and contact guidance from a non-local perspective and the biological relevance of this topic, in this paper we develop a non-local kinetic model describing the combined effect of these two directional cues on the choice of the direction of cell migration. We propose a transport equation implementing a velocity-jump process in which the transition probability describing the choice of the new velocity takes into account the non-local sensing of the fiber network and a chemoattractant. In particular, we propose two possible sensing strategies: in the first case, cells can measure the guidance cues independently and weight them differently; in the second case, cells measure simultaneously the two cues and weight them equally. To the best of our knowledge, this is the first time that a non-local sensing of the fibers distribution defined at a mesoscopic level is considered and is coupled with a non-local sensed chemoattractant. As the kinetic approach allows to derive macroscopic models, we will discuss how the choices of the transition probability, together with the size of the sampling volume and the characteristic variability of the two cues, determine the macroscopic behaviors. This analysis, along with the numerical simulations, will show how different regimes can lead to very different macroscopic equations and that a drift-diffusion model directly stated at the macroscopic level is not always capable of capturing this dynamics.

Specifically, in Sect. 2.1 we present some mathematical preliminaries guiding the reader toward a more clear understanding of the mathematical approach, while the novel characteristic aspects of our transport models are introduced in Sect. 2.2. Then, in Sect. 3 we derive the macroscopic models. Precisely, after the brief introduction of Sect. 3.1 about the well known techniques based on asymptotic expansions we use in the following sections, we shall analyze the possible macroscopic regimes in Sect. 3.2, in relation to both cell size and characteristic lengths of variation of the cues. In Sect. 3.3 and 3.4, we derive the macroscopic limits in the various regimes starting from the transport models implementing the two different sensing strategies. In Sect. 4, we present different numerical tests, in which we numerically integrate the kinetic transport equations in order to qualitatively investigate various scenarios in a two-dimensional setting. We highlight the main features of our approach and make some considerations about the derived macroscopic models. Finally, in Sect. 5, we present an application of our model to an experimental setting, while in Sect. 6 we shall draw some conclusions and illustrate the potential and the possible outcomes of our work.

A Kinetic Model for Chemotaxis on a Fiber Network

Preliminaries

The cell population is described at a mesoscopic level through the distribution density that, for every time and position , gives the statistical distribution of the speeds , where V is the maximal speed a cell can achieve, and of the polarization directions , being the unit sphere boundary in . Thus, the velocity vector will be given by . A macroscopic description for the cell population can be classically (Cercignani 1987) recovered through the definition of moments of the distribution function p. Precisely, we can define the macroscopic cell density in the physical spacethat is the expected mass in . The mass mean velocity that is given by the average of the microscopic velocities of the normalized population can be defined as followsAs a consequence, the first moment of pmay be seen as the mass flow. Concerning higher-order moments, we shall considerSince we are not considering hard spheres, as instead classically done in particle physics, and, in particular, cells are self-propelled entities with an inner energy, we shall not refer to (3) as the energy tensor or momentum flow. In particular, we only consider (3) to be a diffusion tensor of the population that is given by the variance–covariance matrix of the normalized population multiplied by the cell density. The latter gives information about the breadth of the distribution of the microscopic velocities with respect to the average velocity (Painter and Hillen 2018).

The mesoscopic model consists in the following transport equation for the cell distributionwhere the operator denotes the spatial gradient, so that the term takes into account the free particle transport. The term is the turning operator that describes the scattering of the microscopic velocity in direction and speed. In the present case, the typical microscopic dynamics of the cell is the run and tumble (Block et al. 1983; Berg 1983) that prescribes an alternation of runs over straight lines and re-orientations. The choice of the new direction can be random or biased by the presence of external factors, which might attract or repel the cells. The run and tumble is classically modeled by a scattering of the microscopic velocity called velocity jump process (Stroock 1974), characterized by a turning frequency and a transition probability T. The latter models the probability of choosing a certain direction and speed v after a reorientation. For our purposes, we assume that the transition probability does not depend on the pre-reorientation velocity, as classically done in the pioneering work concerning kinetic equations for velocity jump processes (Stroock 1974; Othmer et al. 1988; Hillen 2006) and in the previous works with non-local transition probabilities Loy and Preziosi (2020, 2020). Therefore, the turning operator readsAs T is a conditional probability, it satisfies The mean macroscopic velocity after a tumble is given by the average of Tand the variance–covariance matrixSince we are going to consider two-dimensional bounded domains without loss of cells and no cells coming in, we assume conservation of mass. Thus, we will require that boundary conditions, which must be imposed at the kinetic level, implement no-flux boundary conditions at the macroscopic level (also named biological no-flux boundary conditions) (Plaza 2019). Specifically, we consider, at the kinetic level, a particular prototype kind of no-flux conditions that are specular-reflecting boundary conditions, which means that cells are reflected with an angle of when they hit the wall, as previously done by Loy and Preziosi (2020).

Structure of the Transition Probability

In this section, we introduce the transition probability modeling the decision process of a cell in the presence of two directional guidance cues: a fibrous ECM and a chemoattractant. We consider amoeboid cells moving by contact guidance without proteolysis, , cells hit a fiber and move along the direction of the fiber (Wolf et al. 2003). It has been shown experimentally, for example in the case of glioma cancer cells (Johnson et al. 2009), that randomly disposed fibers imply isotropic cell diffusion, while aligned fibers cause anisotropic diffusion of the cells along the direction of the fibers. The first transport model (the M5 model) for contact guidance was proposed by Hillen (2006), further studied and developed in Painter (2008); Chauviere et al. (2007, 2007), and applied to the study of glioma in Painter and Hillen (2013); Engwer et al. (2015, 2016, 2017); Conte et al. (2020). Following Hillen (2006), we consider a distribution of fibers described as a distribution defined over the space of directions given by the unit sphere in ,that satisfies The last condition means that we are considering a non-polarized network of fibers, so that cells are able to go in both senses in each direction. Being, then, a probability density distribution, we can define the mean direction of the fibersand the variance–covariance matrix of qthat is proved to be the diffusion tensor of the fibers (Hillen 2006). As we consider a non-polarized fiber network, we have thatmeaning that there is no mean direction in the dynamics. When q is a regular probability distribution, the tensor (11) is symmetric and positive definite and, thus, it is diagonalizable. Each eigenvalue represents the diffusivity in the direction of the corresponding eigenvector. This means that, if the eigenvalues are equal, there is isotropic diffusion, while, if they are different, there is a preferential direction of motion, , anisotropic diffusion. This allows to reproduce isotropic/anisotropic diffusion on a non-polarized fiber network (Hillen 2006; Painter 2008). Concerning chemotaxis, we consider a chemoattractant defined in the region by a strictly positive definite functionFor both chemotaxis and contact guidance, we assume that the sensing performed by the cells is non-local, as cells extend their protrusions, through which they sense the environment, up to several cell diameters (Berg and Purcell 1977). The maximum length R of a protrusion is called sensing radius. It has been firstly introduced by Othmer and Hillen (2002) for modeling a non-local gradient of a chemical and, then, used in a number of works for describing the sensing of macroscopic quantities (see Chen et al. 2019 for a review and references therein). In particular, in Loy and Preziosi (2020) and, later, in Loy and Preziosi (2020) it has been used for the non-local sensing of the cues affecting cell polarization and speed.

.

In the present work, we extend the model proposed in Loy and Preziosi (2020) dropping the effect of a cue on cell speed, which will be unbiased, and assuming a double sensing of cues that affect cell polarization. Thus, we assume that both and q are non-locally sensed by a cell that, starting from its position , extends its protrusions in every direction up to the distance R. With respect to a local sensing, assuming a non-local sensing of the fiber network allows to reproduce a wider range of migration strategies that a cell can perform in order to cleverly move toward the chemoattractant. Therefore, we consider the quantitiesNext to the border of the domain , we always consider such that . The term takes into account the possibility of sensing in direction radial fibers starting from the position and up to a distance . For the cell speed, we consider the probability density distribution of the speeds defined on the interval [0, V] that thus satisfies . Precisely, we define its mean speed, energy and variance asrespectively.

Regarding the sensing of the microenvironment, in general cells use different receptors located on their membrane. For instance, the interactions between cells and fibers are mediated by transmembrane receptors, such as the integrins, which allow cell–fiber bindings and the activation of related intracellular pathways. Instead, depending on the chemotactic cue, other receptors or transmembrane units, rather than integrins, can be recruited for the signal transduction. For example, endothelial cells use vascular endothelial growth factor receptors (VEGFR -1 -2 -3) in response to the chemotactic cue given by the vascular endothelial growth factors (VEGFs). Thus, on a first approximation, we can assume that cells respond independently to the information given by these different receptors and the choice of the new direction results from a combination of these two responses. Instead, on a second approximation, we can assume that cells mediate simultaneously the information coming from different receptors and the new direction is the result of a unique weighting process of the measured information. For these reasons, we propose two different transition probabilities to describe these two different sensing strategies: in the first model, the sensing of q and are independent, while in the second model a unique sensing is performed.

In the first model, we introduce a transition probability that is the product of two different independent sensing processesIn this case, a cell located in position extends its protrusions up to a distance R in each direction . Then, the cell measures the field with specific receptors, weighting the related information by , and, independently, it measures the quantity , weighting the related information by . The sensing functions and have compact support in [0, R] and they may be Dirac deltas centered in R, if the cell only measures the guidance cues on its membrane (only on for every ), or Heaviside functions if the cell measures and gives the same weight to q and from to in each direction. Formally, from a stochastic point of view this choice describes an independence of the two phenomena and, thus, the transition probability can be seen as the product of the independent probabilities of q and , , .

The second model prescribes a simultaneous averaging of the guidance cues and q, ,This transition probability describes a cell in position that extends its protrusions and measures in each direction the two quantities and , weighting both simultaneously with the same sensing function . Formally, as the two sensing are not independent, the transition probability cannot be factorized.

In (14) and (15), is a normalization coefficient.

We refer to the transport model (4)-(5) with (14) as non-local independent sensing model, in which the cell averages the two cues independently according to two different sensing functions , . On the other hand, the transport model (4)-(5) with (15) is defined as non-local dependent sensing model, describing cells that sense the two cues at the same time and average them with a unique sensing kernel . In the next section, we briefly present some preliminaries concerning asymptotic limit procedures and, then, we derive the hydrodynamics models in four different scenarios. We compare the macroscopic settings obtained from the two different transition probabilities described in this section.

Derivation of the Hydrodynamic Models

Preliminaries on Asymptotic Limit Procedures

In order to investigate the overall trend of the system, the macroscopic behavior is typically analyzed.

Precisely, we are interested in the evolution of in the emerging regime of the system that may correspond to a diffusion-driven regime or to a drift-driven regime. Therefore, we shall consider a diffusive or a hydrodynamic scaling of the transport equation (4) with (5), respectively, resulting from a proper non-dimensionalization of the system. Diffusive and hydrodynamic limits for transport equations with velocity jump processes have been widely treated by Hillen and Othmer (2000); Othmer and Hillen (2002); Hillen (2006); Loy and Preziosi (2020); Bellomo et al. (2007); Filbet et al. (2005). For the reader’s convenience, in the present section we briefly report a well-known technique based on Hilbert expansions that are aimed at performing the diffusive and hyperbolic limits of our transport models that we shall present in the following sections. Formally, we introduce a small parameter and we re-scale the spatial variable asbeing the macroscopic spatial variable. According to the other characteristic quantities of the system and up to an appropriate nondimensionlization, the macroscopic time scale will becorresponding to the parabolic scaling in the case of a diffusion dominated phenomenon, orcorresponding to the hyperbolic scaling in the case of a drift driven phenomenon. As widely treated by Othmer and Hillen (2002), we consider a framework in which, up to the spatial scaling (16), we can expand the transition probability asThis means that there are different orders of bias. If we assume that , we denote by and the corresponding operators defined by and , respectively, and we assume that and The corresponding first and second order moments are given byandConsidering the Hilbert expansion of the distribution function pif there is conservation of mass, we have that all the mass is in (Hillen and Othmer 2000), ,where . Furthermore, for performing the diffusive limit we shall assume that (Hillen and Othmer 2000). The fundamental property for performing the diffusive limit requiresmeaning that the leading order of the drift velocity vanishes. This is coherent with the fact that the time scale is chosen because macroscopically the phenomenon is diffusion-driven.

Equation (4), rescaled according to (16)-(17), readswhere represents the turning operator defined by . Equating the terms of equal order in , we obtain the following system of equations.

In :In :In :Eq. (26) impliesthat is the equilibrium state of order zero. From Eq. (27), by inverting , we getPrecisely, the functional solvability condition necessary for inverting iswhich is satisfied because (24) and (19b) are satisfied.

Integrating (28) over , we obtain the macroscopic diffusive limit given by (dropping the dependencies)being the diffusion motility tensor. Equation (32) is a diffusion–advection equation, where is the drift velocity of first order. If (24) does not hold, a hyperbolic scaling is more appropriate. In this case, rescaling the variables as in (16) and (18), the transport equation (4) readsSince the equilibrium state is the same as before, we consider the Chapman–Enskog expansion of p in the formwhere . Substituting (34) in (33) and integrating the equation at the order over , we obtainThis is an advection equation modeling a drift driven phenomenon. We address the reader to Loy and Preziosi (2020) for further details about the derivation of the macroscopic model. Concerning the boundary conditions (Plaza 2019), at the macroscopic level no flux boundary conditions for the diffusive limit readbeing n the outward normal to the boundary, whilst for the hyperbolic limit they are given by

Mesoscopic Analysis of Two Non-Local External Cues

In line with the analysis proposed by Loy and Preziosi (2020, 2020), in order to qualitatively evaluate the effects of the non-locality at the macroscopic level, we study the impact of both directional cues and q with respect to the size of the cell (which is related to the cell sensing radius R). Precisely, the following analysis refers to spatially heterogeneous fiber networks and chemicals, as for homogeneous cues there is a loss of significance in considering a non-local sensing. As in Loy and Preziosi (2020), we introduce the characteristic length of variation of It allows to approximate with a positive quantitywhere we neglect higher-order terms in . Beside the above defined characteristic length of variation of the chemoattractant, in the present work we define an analogue quantity for the fibers distribution. We define in particularIn this case, we can approximate with a positive quantityThis definition of takes into account the variation in directionality of the fibers in space, which is what actually influences cell orientation, more than the spatial variation in the density of the ECM. We analyze the possible scenarios depending on the relation between R, and . For this purpose, let us introduce the parametersandthat quantify the cell measuring capability with respect to the characteristic lengths of variation of the guidance cues and q. In particular, , means that the sensing radius is smaller than the characteristic length of variation of q (, respectively) and, thus, a single instantaneous sensing of the cell is not capable of catching the total spatial variability of q (, respectively). On the other hand, if , the sensing radius is large enough to capture the spatial variability of q (, respectively). If we consider the two cues separately, in the first case we expect that the sensing of q (, respectively) induces a diffusive behavior, while in the second scenario the overall behavior induced by q (, respectively) is drift driven. However, as we are considering the two guidance cues simultaneously affecting cell polarization, we have to take into account four limit cases: We remark that the concept of fast and slow variation refers to the cell sensing capability and it is intended with respect to the relation between cell sensing radius and characteristic variation length of the cues.

fast variation of both external cues: ;

slow variation of both external cues: ;

fast q and slow variation: ;

fast and slow q variation: .

In the next section, we analyze the macroscopic trends in the cases . Precisely, in case (i), a Taylor expansion for both cues cannot be used, since there is no guarantee that the first order approximations are positive, as well as in case (iii) and (iv) for q and , respectively. A priori, in order to quantify the relative contribution of chemotaxis to contact guidance, we introduce the parameterthat is larger than 1 if contact guidance prevails, whilst it is smaller then 1 if chemotaxis is stronger. Due to (41) and (40), we have that, despite its definition, does not depend on the size and sensing capability of the cell, but only on the characteristics of the cues, as . In particular, if is larger than , , , it means that the gradient of q is steeper than the one of , thus enhancing a stronger effect of contact guidance on the dynamics. We may also observe that in the case of fast q and slow variation (case (iii)) we always have while in the case of fast and slow q variation (case (iv)) we always have , , contact guidance is weaker than chemotaxis.

Amoeboid Motion and Chemotaxis: Non-Local Independent Sensing

We first consider the non-local independent sensing case (4)-(5) with (14). We recall the expression of the transition probabilityThe average of T, which is the equilibrium velocity of the cell population, is given by

Case (i): Fast Variation of Both External Cues

In this case, we shall chooseAs a consequence of the fact that T cannot be expanded in powers of after the re-scaling with (16), we have that the equilibrium velocity is directly given by (43). Therefore, we perform a hyperbolic scaling that leads to the following macroscopic equation for the cell density:with given by the re-scaling of (43) with (16).

Case (ii): Slow Variation of Both External Cues

In this case, we can expand both and and consider the approximations (37) and (39) for . Thus, we approximate the transition probability by substituting (37) and (39) in (14), and we obtain the following approximation for the turning kernel where we neglect higher orders terms in . In the latter,andThe quantities are the weighted (by ) measures of the sensed linear tracts in every direction, whilst are the averages of on [0, R]. Then, we introduce the small parameterand re-scale the space variable as , gettingThis expression means that the equilibrium is determined by the fiber distribution. Moreover,whereIt can be easily checked that and verify (19a) and (19b), respectively. Because of (12) and (47), we get , meaning that we are in a diffusive regime. The diffusive limit leads to the advection-diffusion equation (32). The explicit form for the zeroth-order macroscopic diffusion tensor iswhile for the macroscopic first-order velocity isTherefore, the resulting diffusion–advection equation reads (dropping the dependencies)whereare the sensitivities. The diffusion represented by the cell motility tensor (48) only depends on the fibers distribution, while the advective term has two contributions differently weighted by the sensitivities in (51). We remark that, in the diffusive regime we obtain the same macroscopic behavior postulated by Keller and Segel (1970), with the logarithmic chemotactic sensitivity given in (51). The term depends on both the fibers distribution and the chemotactic field; it never vanishes if is not the null vector, since it may be proved that is invertible. In the isotropic case, corresponding to randomly disposed fibers, , when is proportional to the identity matrix, then is parallel to , which, thus, represents the anisotropy direction. On the other hand, when is anisotropic and is not parallel to the leading eigenvector of , then the migration does not follow the dominant direction of the fibers, but rather its projection on . The second contribution in the drift term, , , is a measure of the velocity field induced by the spatial variation of the fiber distribution, which determines the cell microscopic velocities. This term vanishes if the fibers distribution is homogeneous in space. However, if q is homogeneous in space, even in case of competing cues, , , in general the advective term does not vanish. Instead, in case of cooperating cues, , when is an eigenvector of with eigenvalue , migration is in direction with a kinetic factor . In intermediate scenarios, migration happens in the projection , but, if q is not homogeneous, the dynamics is more complex and, even in case of cooperation, we cannot conclude anything about additivity effects.

Case (iii): Fastqand SlowVariation

In this case, we can only expand with Taylor series the chemoattractant, as in (37), and the turning kernel (14) can be approximated aswhereandare both different from zero. In this case we may chooseand, re-scaling the space variable as (16), we get . Let us now introducea new probability density distribution describing a non-local average of the fibers distribution according to the sensing kernel and normalized by the measure of the sensed linear tract over the direction . With this notation, we can definewhich is the corresponding diffusion tensor of the fibers, andwhich is their corresponding mean direction. Precisely, both are defined taking into account the whole neighborhood sensed by the cells. In this case, does not vanish in general in , as it is given byTherefore, we perform a hyperbolic limit that leads to the drift equation (35) with mean velocity (56). Precisely, this mean velocity is related to a non-local average of the diffusion tensor of the fibers projected on , and a non-local average of the mean fiber direction depending on the local chemoattractant . In this case, there is not an evident additivity effect of the two cues and several scenarios are possible.

Remark In the particular case in which also the new introduced probability density has the same symmetry property Q3 (described in Sect. 2.2), or if we consider (local sensing of the fibers), a parabolic scaling could be performed instead of a hyperbolic one. This would lead to a macroscopic diffusion–advection equation with mean velocity . Without chemotaxis, we would recover the classical model for contact guidance (Hillen 2006), which gives rise to a fully anisotropic diffusive equation at the macroscopic level. The presence of a non-local chemoattractant, even when , is responsible for the emergence of the drift correction term.

Case (iv): Fastand SlowqVariation

This last case only allows for the Taylor expansion of the distribution function q, as in (39). Therefore, the turning kernel can be approximated aswhereandboth different from zero. In this case we may chooseand, by re-scaling (57) with (16), we get . Hence does not vanish in , as it is given byand the macroscopic equation is given by (35). The mean velocity (58) is a linear combination of a non-local measure of the chemoattractant over the fiber network and a non-local measure of weighted by the directional average of the spatial variability of the fiber direction.

Remark

If we consider a local sensing for the chemoattractant, , , we obtain a macroscopic advection-diffusion equation. Here the macroscopic velocity would be induced by the spatial variation of the fiber direction distribution , and the measure of does not affect the choice of the direction. In this case, if vanishes, the model would reduce to a fully anisotropic diffusive equation (Hillen 2006).

Amoeboid Motion and Chemotaxis: Non-Local Dependent Sensing

Concerning the non-local dependent sensing case (4)–(5) with (15), we recall the expression of the transition probabilitywithThe macroscopic velocity is here given byThe macroscopic limits can be performed analogously to the previous section and using the techniques described in Sect. 3.1. The choice of the parameter for the cases will be the same of the non-local independent sensing model, since it does not depend on the kind of model (independent or dependent sensing), but only on and , i.e., on the rapidity of variation of and q with respect to the cell sensing capability.

In this case, we cannot consider the expansions (39) and (37), and, thus, we cannot expand the turning kernel. Its non-vanishing average is given by (59). Therefore, we perform a hyperbolic limit leading to (35) with macroscopic velocity (59).

When, instead, the maximum sensing radius R is smaller than both characteristic variation lengths, we consider the positive expansions (39) and (37) and substitute them into (15). Neglecting the higher order terms in , we get the approximationwithandRe-scaling the space variable as in (16), we findandwithIt can be easily checked that and verify (19a) and (19b), respectively. Therefore, , because of (12), and we can perform a diffusive scaling that leads to the zero-order macroscopic diffusion tensorand to the macroscopic first-order velocityThe macroscopic advection-diffusion equation (32) now reads (dropping the dependencies)whereSimilar observations to the case (ii) of the non-local independent sensing model may be done, except that, here, there is a unique sensitivity that weights equally the two contributions in the advection term (62).

In this case, we expand only the chemoattractant , as in (37), and the turning kernel (15) can be approximated aswithandboth different from zero. Re-scaling the space variable as (16), we find . Let us defineandmean fiber direction and diffusion tensor, respectively. They are related to the average, through the function , of the non-local fiber distribution evaluated on the whole neighborhood sensed by the cells. Hence, in this case, does not vanish in general in , as it is given byTherefore, the macroscopic advection equation has an expression analogous to (35) with macroscopic velocity (67), which represents a linear combination of a non-local weighted average of the diffusion tensor of the fibers projected on , and a non-local average of the mean fiber direction depending on the local chemoattractant . As for the independent sensing case, there is not a simple additivity effect of the two cues and several scenarios are possible.

In this case, again, we can only consider the positive approximation (39), and the transition probability rewrites aswhereandboth different from zero. As before, by re-scaling (68) with (16), we get and we have that the average velocity . In particular, it is given byand, thus, we perform a hyperbolic limit leading to (35). The mean velocity (69) is a linear combination of a non-local measure of the chemoattractant over the fiber network and a non-local average of weighted by the directional average of the fiber direction spatial variability.

Comments

We can observe that, if , the two non-local transport models for independent and dependent sensing are equal, while, if the sensing kernels are not Dirac deltas (even if ), the transport models are always different. Instead, at the macroscopic level, with any choice of the sensing functions the models coincide only in case (ii), i.e., when we have a slow variation of both cues. In this case, in fact, the macroscopic limits are different only if , while in the cases (iii), i.e., when we have fast q and slow variation and iv), i.e., when we have fast and slow q variation, they are different if the sensing kernel are not Dirac deltas (even if ). The relevant difference concerns the macroscopic transport velocities (see (56) and (67) for the case (iii) and (58) and (69) for the case (iv). In fact, in the cases (iii) and (iv), for the non-local dependent sensing model, since only one cue is considered non-locally and both cues are averaged with the same sensing function , we have a weighted average on of the non-local quantities, which results in the weighted averages in the second terms of (67) and (69). Table 1 presents a general summary of all the models we derived, while we summarized these remarks in Table 2.

Table 1

Summary of the models (dropping the local dependencies in )

Case	Non-local independent sensing (4)-(5)-(14)	Non-local dependent sensing (4)-(5)-(15)
i) fast S and slow q	Drift dominated	Drift dominated
	U_T=cU¯∫_S^d-1v^∫0Rγ_S(λ)S(ξ+λv^)dλ∫0Rγ_q(λ)q(ξ+λv^,v^)dλdv^	U_T=cU¯∫_S^d-1v^∫0Rγ(λ)S(ξ+λv^)q(ξ+λv^,v^)dλdv^
ii) slow S and q	Drift-diffusion	Drift-diffusion
	DT0=DD_q	DT0=DD_q
	UT1=U¯Γ^q∇·D_q+Γ^SD_q∇SS	UT1=U¯Γ∇·D_q+D_q∇SS
iii) fast q and slow S	Drift dominated	Drift dominated
	U_T=U¯Γ0Sc_0S∫_R_+γ_q(λ)∫_S^d-1v^q(ξ+λv^,v^)dv^dλ	U_T:=U¯c_0∫_R_+γ(λ)∫_S^d-1v^q(ξ+λv^,v^)dv^dλ
	+U¯Γ1Sc_1∇S∫_R_+γ_q(λ)∫_S^d-1v^⊗v^q(ξ+λv^,v^)dv^dλ	+U¯c_1∇S∫_R_+λγ(λ)∫_S^d-1q(ξ+λv^,v^)v^⊗v^dv^dλ
iv) fast S and slow q	Drift dominated	Drift dominated
	U_T=U¯Γ0qc_0∫_S^d-1v^q∫0Rγ_S(λ)S(ξ+λv^)dλdv^	U_T:=U¯c_0∫_S^d-1v^q∫0Rγ(λ)S(ξ+λv^)dλdv^
	+U¯Γ1qc_1∫_S^d-1v^⊗v^∇q∫0Rγ_S(λ)S(ξ+λv^)dλdv^	+U¯c_1∫_S^d-1v^⊗v^∇q∫0Rλγ(λ)S(ξ+λv^)dλdv^

Table 2

Summary of the comparison of the models for different choices of the sensing functions. indicates the cases in which the models coincide, while the ones in which the models are different

	γ_q=γ_S=γ=δ	γ_q=γ_S=γ≠δ	γ_q≠γ_S
Meso models (4)-(5)-(14) and (4)-(5)-(15)	✓
Macro models case (i): fast S and q	✓
Macro models case (ii): slow S and q	✓	✓
Macro models case (iii): fast q and slow S	✓
Macro models case (iv): fast S and slow q	✓

Numerical Simulations

We now present simulations of the kinetic transport model (4)-(5) for non-local independent sensing (14) and non-local dependent sensing (15) in order to show some key model features. Precisely, we numerically integrate the transport equation to approximate the density distribution p, as in Filbet and Yang (2014); Loy and Preziosi (2020), and in turn the corresponding macroscopic density via (1). For computational convenience, we restrict to a rectangular 2D region .

We specify q using the standard circular distribution given by the bimodal von Mises distribution (Mardia and Jupp 2009)where is the modified Bessel function of first kind of order andIt can be proved that the first trigonometric moment is (Hillen et al. 2017), and, therefore, is the mean direction of the fibers located at point (Mardia and Jupp 2009). This function also satisfies Q3 and its variance–covariance matrix is given by Hillen et al. (2017)where is the identity tensor in , while k and are functions of . Moreover, the circular-variance is given by the scalarthat represents the degree of alignment of the fibers at point .

We propose three sets of numerical tests in which we integrate the kinetic transport equation (4) and visualize the macroscopic density (1):

a preliminary test for showing the effects of a non-local chemoattractant in the presence of a local fiber network;

a comparison between the two sensing strategies with different sensing kernels;

a set of simulations in different scenarios (, for different ), allowing to make some observations about the differences between the emerging macroscopic regimes.

Test 1: Local ECM Sensing and Non-Local Chemotaxis

Test 1 is designed to highlight the effect of the presence of a chemoattractant on the behavior of cells migrating by contact guidance over a non-polarized fiber network and evaluating locally its alignment. This means that for this test we consider . In the absence of additional cues, when fibers are sensed locally, even though the distribution q is anisotropic, cells are distributed symmetrically in each given direction as there is no drift induced by an asymmetric external cue (Hillen 2006; Painter 2008). Conversely, in presence of a chemoattractant, if it is sensed non-locally, a chemotactic velocity is imposed. Formally, we are dealing with (14) in which . In particular, we consider a regionwith and in which the fibers are strongly aligned along the direction identified by . In particular, for , , such that . In the rest of the domain fibers are uniformly distributed. The chemoattractant has a fixed Gaussian profileIn particular, for Test 1 we chose . The initial condition for the cell population is a Gaussianwith and . In this first test, the initial condition for the cell population is centered in , , the center of the region (see Fig. 1a). We remark that without chemoattractant, cells would diffuse anisotropically in the preferential direction of the fibers , forming the well known ellipsis (Painter 2008), which represents cells moving with the same probability along directions and . In the present case, because of the presence of a chemoattractant, the symmetry is broken, and, even if q describes a non-polarized fiber network, there is a preferential sense of motion in direction (see Fig. 1d–f). Precisely, cells migrate along the fibers in the direction identified by , corresponding to the preferential sense imposed by the presence of the chemoattractant located in the upper-right corner of the domain . Given this directional setting, the cell population dynamics is also greatly affected by the strength of the chemoattractant, which depends on and , by the degree of the alignment , which depends on k(x, y), and by the sensing radius R.

Fig. 1

Test 1: evolution of the initial cell distribution in (a) for local q and non-local with sensing function . The sensing radius of the cells is set to , while is (72) with and . (b): initial cell orientations. (c): trajectory of the center of mass of the cell population, where each black dot is plotted every . (d–f): evolution of the macroscopic density. (g–i): polarizations of the cells (color figure online)

Another important aspect is the sensing function , which influences the transient dynamics and, especially, the relaxation time. For example, in the case of a Heaviside function, the relaxation time is twice the relaxation time needed when is a Dirac delta (see also Loy and Preziosi 2020). Moreover, we analyzed the average cell polarization at every position , given by the mass flow (2), that is the average velocity times the macroscopic cell density. The cell microscopic orientations are initially randomly distributed and they start from a vanishing initial speed (see Fig. 1b). Then, they start to orient and align along the fibers, migrating upward in the direction individuated by the angle , since cells sense the chemoattractant (see Fig. 1g–h). Eventually, when cells reach the level , the microscopic directions polarize toward the chemoattractant (see Fig. 1i). The center of mass plotted in Fig. 1c stays in the region during cell migration along the fibers bundle in , and it moves out of only when it reaches . The black dots are plotted every and it is clear that the highest acceleration happens when cells are on the bundle of fibers, while they are slowed down when they start to move out of the fibers stripe .

Test 2: Non-Local ECM Sensing and Chemotaxis

In this second test, we compare the non-local independent sensing model and the non-local dependent sensing model. We assume fibers to be distributed similarly to the previous test, , fibers are highly aligned in , defined this time by and (see Fig. 2b). Here, for , , that corresponds to , and . In the region fibers are uniformly distributed. The initial condition of the cell population is (73) with (see Fig. 2a) while the chemoattractant has a fixed profile located as in Test 1, with and . We compare the dynamics of the cell population in the following four settings: Results of these simulations are shown in Fig. 2.

Fig. 2

Test 2 Time evolution of the initial distribution given in Fig. 2a in the four settings 1-4. The sensing radius of the cells is and is (72) with and . Setting 1 is represented in Figs. (c–f): local q and non-local , . Setting 2 is represented in Figs. (g–j): non-local q and with sensing functions . Setting 3 is represented in Figs. (k–n): non-local q and , independent sensing with . Setting 4 is represented in Figs. (o–r): non-local q and , dependent sensing with (color figure online)

local fiber distribution and non-local chemoattractant, as in Test 1, , (14) with and ;

non-local sensing with Dirac Deltas for both q and , , we integrate (4) with (14) or (15) with ;

non-local independent sensing with Heaviside sensing functions for both and q, , (4)-(14) with ;

non-local dependent sensing with Heaviside sensing function for q and , , (4)-(15) with .

We can observe that, in all settings, cells start from (1, 0.5), they are attracted by the chemoattractant and, on their way towards , they cross the aligned fibers region , climbing up this region in the direction . Eventually, in all the cases, cells reach the chemoattractant, but the dynamics, as well as the transient time, is influenced by different sensing kernels and the local or non-local sensing strategy, even though the differences are not extremely evident. Although settings 3 and 4 in Fig. 2, which are related to the case of independent and dependent cues, respectively, do not seem to show very strong differences, we can observe that in case 3 (see Fig. 2k–n) the tendency of going in both the directions , determined by q, and , determined by , appears more marked because of the independent sensing. In contrast, this behavior is less evident in case 4 and it results the least evident in the case in which cells deal with a local sensing of the fibers (setting 1), with also a general slow down of the dynamics.

Test 3: Non-Local Independent Sensing Model for the Comparison of Cases

Test 3 is designed to explore the extent to which the macroscopic cell behavior changes depending on the relation between R, and . Precisely, we compare cases (i), (ii), (iii) and iv) according to different mutual variations of the external cues q and . The aim here is to illustrate the importance of choosing the most appropriate macroscopic model, as the aggregate cell behavior can be very different according to the regime prescribed by the system parameters. We perform the analysis for the non-local independent sensing case with , as it is a case in which the transport model is different from the dependent sensing model. Moreover, the independence of the two sensing allows to visualize more efficiently the two distinct directional effects (contact guidance and chemotaxis), as observed from the results of Fig. 2.

We consider the turning kernel describing contact guidance with q given by (70), , and coefficient k(x, y), modulating the strength of the alignment, given by the gaussian distributionHere, and (Fig. 3d). This choice describes a setting in which fibers are more aligned in the central circular region of the domain and uniformly disposed in the rest of it. We consider different values of in order to obtain different values of : corresponds to and corresponds to . Details about the estimation of for a Bimodal Von Mises distribution of fibers q are given in Appendix A. The chemoattractant has the fixed profile (72) with and . In the simulations, we consider three different values for the chemoattractant variance in order to obtain different values of : , which corresponds to in Fig. 3a; , which corresponds to in Fig. 3b; , which corresponds to in Fig. 3c. The initial cell distribution for all the settings presented in Figs. 4-8 is given by (73) with , , . Precisely, we present five sets of simulations that are summarized in Table 3.

Fig. 3

Test 3 Three different chemoattractant distributions used for comparing models . The chemoattractant profile is given by (72) with and (a) , corresponding to , (b) , corresponding to , and (c) , corresponding to . The fibers distribution in sketched in (d) (color figure online)

Fig. 4

Test 3 Case (i) (fast variation of both cues) with non-local q and , sensed with an independent sensing through the kernels . is given in Fig. 3a with and , so that . The fibers distribution q has a space dependent parameter k given by (74) with , so that . The sensing radius of the cells is (color figure online)

Fig. 8

Test 3 Case (iv) (fast and slow q variation) with non-local q and , independent sensing with . is given in Fig. 3a that corresponds to , while for the fiber distribution , so that . The sensing radius of the cells is (color figure online)

Table 3

Summary of the simulations presented in Test 3

l_S	l_q	R	Case	η	Figs.
0.002	0.0031	0.7	(i)	<1	4
0.25	0.0031	0.7	(i)	≫1	5
0.055	0.031	0.02	(ii)	>1	6
0.25	0.0031	0.2	(iii)	≫1	7
0.002	0.031	0.02	(iv)	<1	8

In Fig. 4, we consider the case in which , , fast variation of both external cues (case i). This corresponds to a strongly hyperbolic macroscopic behavior with macroscopic velocity given by (43). In Fig. 4 we can observe that cell behavior is not diffusive and the cluster of cells is quite compact. Moreover, when cells reach the central region where fibers are strongly aligned in the direction (as shown in Fig. 3d), which is perpendicular to the favorable direction induced by , they surround this region and go over it toward . In this setting, the parameter defined in (42) is slightly smaller then 1 and, in fact, chemotaxis prevails in the overall dynamics, as the stationary state is clearly peaked on the chemoattractant profile, but the fibers structure influences the transient.

In Fig. 5, we consider with and, consequently, (see Fig. 3c). Concerning the fibers, we have , so that , and the sensing radius is . This setting falls again in case (i), fast variation of both cues, but the behavior is different with respect to the previous simulation in Fig. 4. The chemoattractant in Fig. 3c, in fact, is spread over the whole domain and, actually, the quantity is almost times the considered in Fig. 3a (and used for the simulation in Fig. 4). Even though we are still in a hyperbolic case and cells are guided by the strong drift (43), as R is slightly larger than and is large, the cell cluster diffuses a bit more in the domain. When cells reach the region of strongly aligned fibers, they start to surround it (see Fig. 5a–c), but, as , some of them do not surround the region, slow down and partially tend to align along the fibers. In Fig. 5c, for instance, we have a high density of cells in the strongly aligned fiber region. Eventually, cells manage to overcome the area of highly aligned fibers and they tend to converge to the chemoattractant profile (see Fig. 5d). In this setting, the overall dynamics is greatly affected by the fibers and, in fact, .

Fig. 5

Test 3 Case (i) (fast variation of both cues) with non-local q and , independent and sensing with . is given in Fig. 3c that corresponds to , while for the fiber distribution , so that . The sensing radius of the cells is (color figure online)

The second scenario, illustrated in Fig. 6, refers to case (ii), slow variation of both cues, since the sensing radius is smaller than both and . At the macroscopic level, the behavior of the system is described by the diffusion–advection equation (50) with macroscopic velocity (49). Actually, in Fig. 6 we can observe a highly diffusive behavior, as the cell macroscopic density has invaded almost half of the domain before even starting to feel the influence of the fibers. If we compare the same time step in Figs. 5b and 6b, we see how cells are reaching in both cases the fibers and feeling the region in which fibers are aligned the most. However, in Fig. 5b the cell cluster is much more compact than in Fig. 6b, where, instead, cells already occupied half of the domain, because of diffusion, and we have a high density of cells both in the region that is close to the strongly aligned fiber region and around the initial position. Then, cells start surrounding the central region of strongly aligned fibers, because they already sense the chemoattractant, and, once they have overcome this area, they tend to the chemoattractant profile (see Fig. 6b–d). In particular, in the transient time, cells accumulate the most at the sides of the region with highly aligned fibers. In this specific setting, and, in fact, contact guidance highly affects the dynamics.

Fig. 6

Test 3 Case (ii) (slow variation of both cues) with non-local q and , independent and sensing with . is given in Fig. 3b that corresponds to , while for the fiber distribution , so that . The sensing radius of the cells is (color figure online)

The third scenario, illustrated in Fig. 7, refers to the case iii) describing fast q and slow variation, since the sensing radius is smaller than , but it is larger than . The macroscopic setting is described by a drift dominated equation with drift velocity given by (56). As , we have that the chemoattractant induces a strong diffusivity, but being , the alignment of fibers strongly affects the dynamics (see Fig. 7c). Comparing, in addition, Figs. 6b and 7b, we have now that the highest cell concentration is in the mean fiber direction in the region surrounding the center of the domain, where the fibers are aligned with a higher degree. As already observed in Sect. 3, this scenario prescribes and, in fact, contact guidance dominates again the dynamics.

Fig. 7

Test 3 Case (iii) (fast q and slow variation) with non-local q and , independent and with sensing function . is given in Fig. 3c, so that , while for the fiber distribution , corresponding to . The sensing radius of the cells is set to (color figure online)

Eventually, for a sensing radius smaller than , but larger than , we are in the case of fast and slow q variation and the macroscopic behavior is approximated by a hyperbolic equation with drift velocity given in (58). Results of the corresponding simulation are presented in Fig. 8. Here, the chemoattractant has the fixed profile shown in Fig. 3a. Cells diffuse in the domain because is smaller than 1, and they start moving in a region with randomly disposed fibers (see Fig. 8a). Then, they mainly follow the preferential direction thanks to the presence of the chemoattractant. In fact, it induces a strong drift because of the high non-locality, determining . Here chemotaxis is slightly dominating the dynamics and, in fact, .

Lastly, we remark that, due to different parameter settings, the relaxation time toward the final equilibrium configuration is different in the scenarios illustrated in Figs. 4-8. Thus, in each figure, we reported the same intermediate time steps for all the presented scenarios (corresponding to the subfigures (a)–(c)) and a different final time step (corresponding to the subfigure (d)) in order to catch the final equilibrium configuration. Precisely, we observe that in Figs. 4 and 8 the time variation is much faster with respect to the other cases. In these scenarios, in fact, we have , i.e., chemotaxis is prevailing and cells move faster toward the equilibrium configuration. Instead, in Figs.  5-7 we have , which means a stronger influence of the contact guidance phenomenon. In particular, the fiber orientation is competing with the preferential direction of movement given by the chemoattractant and cells dynamics are slowed down, leading to much larger relaxation times. For better clarity, videos of the simulations of Test 3 are included as Supplementary material.

Application: Cell Motility on a System of Electrospun Fibers Under a Gradient of VEGFs

As it was mentioned in Introduction, an important aspect related to the study of multi-cue environments concerns the design of engineered scaffolds allowing to direct cell migration by multiple directional cues. In fact, understanding the controlling factors of cell migration and the interplay between them is necessary for designing implants with optimal cellular infiltration and implant integration with native tissue. Moreover, it is important for performing experiments related to phenomena that occur in vivo in physiological and pathological processes. In this section, we aim at showing how our theoretical framework for the description of a double cue environment guiding cell orientation is actually able to replicate the experimental results presented by Sundararaghavan et al. (2013). In this work, the authors consider a model system of electrospun hyaluronic acid fibers, which is an engineered scaffold that is able to mimic native tissue and have the control over fiber orientation. On this scaffold, they evaluate the motility of human umbilical vein endothelial cells (HUVECs) considering the cases of aligned and non-aligned fibers of hyaluronic acid in the presence of a gradient of vascular endothelial growth factors (VEGFs). This setting reproduces a quite common biological situation. Endothelial cell migration towards increasing gradients of VEGFs, in fact, is a well-known characteristics of the process of tumor angiogenesis. Especially in hypoxic situation (i.e., when there is a lack of oxygen in the tumor microenvironment), tumor cells are known to produce and release growth factors in the environment in order to attract endothelial cells and start the process of blood vessel formation or remodeling (Onishi et al. 2011). When the underlying tissue on which endothelial cells are migrating is characterized by aligned fibers, the chemotactic cue given by the VEGF gradient has to be integrated with the guidance response to the fiber network. This is, for instance, the case of brain tumors, where the endothelial cell response to tumor-produced VEGF gradients is influenced by the strong alignment of the fiber structure characterizing the brain tissue (e.g., see Lamalice et al. (2007) and reference therein).

In this section, we want to show how our model is able to recover different behaviors of the cells in response to two main experimental settings: fibers parallelly aligned to the direction of the chemotactic gradient or fibers perpendicularly aligned to the chemotactic gradient. Precisely, we would like to capture the cell dynamics shown in Fig. 3C and D of Sundararaghavan et al. (2013). As initial distribution of cells, we consider (73) with , , and , while the chemoattractant is given by the linear functionwith . The fibers are described with (70) and they are strongly aligned along the direction identified by in the case represented in Fig. 9a, while they are strongly aligned along the direction given by in the case in Fig. 9b. In particular, in both cases, , so that . Moreover, in both scenarios we quantify the cell mean speed along the direction asHere, represents the physical space computational domain and is the weighted -space with weight function .

Fig. 9

Application: migration of cells on a system of electrospun fibers under a VEGF gradient. Cell macroscopic density at time in response to a chemotactic gradient oriented toward the north (as indicated by the arrow) is illustrated in two scenarios: fibers horizontally oriented in the entire domain (a) and fibers vertically oriented in the entire domain (b). The cell sensing radius in both cases is . The time evolution of the cell mean speed defined in (75) is shown in (c) for the case of perpendicular cues (continuous line that corresponds to scenario (a)) and parallel cues (dashed line that corresponds to scenario (b))(color figure online)

Figure 9 shows the resulting cell behavior in the two scenarios. In particular, in Fig. 9a, we consider the case of perpendicular cues, meaning that the angle between the mean fiber orientation and the chemotactic gradient direction is . In this case, the competition between the two cues is evident. In fact, due to the strength of the fiber alignment, cells starting from their initial gaussian distribution align along the fiber direction without being able to move toward the chemoattractant gradient. This is in line with the results shown in Fig. 3C of Sundararaghavan et al. (2013), where cells are seen to move along the direction of fiber alignment with apparently no influence of the VEGF gradient. Looking at the mean speed evolution in Fig. 9c (continuous line), we observe that the cell mean speed along direction rapidly decreases toward zero, meaning that the migration of cells in the direction of the increasing chemotactic gradient is hampered. In Fig. 9b, instead, we analyze the case of parallel cues, i.e., the angle between the mean fiber orientation and the chemotactic gradient direction is 0. The cooperation between the two cues here is evident. In fact cells rapidly migrate following the aligned fibers toward the regions of greater chemoattractant concentration. The numerical results are in a good agreement with the experimental ones, shown in Fig. 3D of Sundararaghavan et al. (2013), where a bias of cells moving up the gradient and along the fiber orientation is observed. Concerning the time evolution of the cell mean speed along direction (dashed line in Fig. 9c), there is a fast increase of as long as cells are moving toward the upper part of the domain where there is the highest computation of the chemoattractant. Once cells have reached the upper border, the mean speed does not increase anymore. In fact, as an effect of the boundary conditions implemented in the simulations and the non-locality (through the cell sensing radius), cells will stop and accumulate in the region of greater chemoattractant.

Conclusion

In this work, we have proposed a kinetic model for describing cell migration in a multi-cue environment. Moreover, we have considered that cells perform a non-local sensing of the environment up to a distance R (named the sensing radius) from its nucleus. Concerning the environmental stimuli, in the present model there are two guidance cues affecting cell polarization, and, thus, cell direction of motion: contact guidance, that is, a bi-directional cue, and a chemical gradient, that is, a mono-directional cue. We remark that, to the best of our knowledge, this is the first time that a non-local sensing in the physical space of the mesoscopic distribution of fibers is considered.

We introduced two novel classes of models: in the first one, cells perform an independent sensing of the fibers and of the chemical in their neighborhood, while, in the second one, cells average the chemical and the fibers with the same sensing kernel. In the two cases, a particular attention was devoted to the identification of the proper macroscopic limit according to the properties of the turning operator. In presence of a heterogeneous network of fibers and chemical, we detected two parameters, and , that measure the relation between the cell sensing radius (R) and the characteristic lengths of variation ( and ) of the two cues. These parameters allow to discriminate between a diffusion-driven regime with an advective correction and a drift-driven regime. Precisely, when the sensing radius does not exceed the characteristic length of the cues, the bi-directional nature of the fiber network prescribes a diffusive regime, otherwise the hyperbolic scaling leads to a macroscopic drift driven regime. We also defined a new parameter that is independent on cell size or sensing capabilities and quantifies the relative contribution of contact guidance to chemotaxis. It provides a first separation between the cases of fiber-dominating and chemotaxis-dominating dynamics ( or , respectively). A common feature we noticed in different cases is the resulting dependency of the macroscopic velocity on both the fiber network and the chemoattractant. This aspect enhances the non-trivial influence of contact guidance on the cell drift and this interdependence is in accordance with the model proposed by Wagle and Tranquillo (2000). Moreover, in absence of the chemoattractant, the fiber impact on the drift term could persist for spatial heterogeneous fiber distributions, in accordance to what is observed by Hillen (2006). This feature represents a step forward with respect to Wagle and Tranquillo (2000), in which the drift is a function of contact guidance only through to the presence of a chemical gradient (meaning that, without chemoattractant, there would be no drift).

The numerical simulations of the transport equations pointed out the main features characterizing the two classes of models and the possible scenarios that they are able to capture. We observed that the presence of two cues influencing cell polarization, even when the fibers are sensed locally (Test 1), ensures a preferential sense of motion for cells laying on regions of highly aligned non-oriented fibers, and the non-locality enhances this behavior (Test 2). Moreover, these non-local aspects brings a further level of detail to our model, allowing to obtain different macroscopic behaviors depending on the characteristics of the two sensing (Test 3). We did not observe remarkable differences between the independent and the dependent sensing models, when we assume in the former the same sensing kernel for fibers and chemoattractant (, when ). However, if there are biological observations sustaining the possibility that a cell might implement different strategies for sensing the fibers and the chemoattractant, it would be possible to use our framework (in its independent sensing version) to investigate this scenario. This could also allow to compare the possible outcomes of different sensing approaches with the case of a unique and common sensing strategy. Moreover, Test 3 showed the importance of deriving macroscopic equations from the underlying microscopic dynamics and in the appropriate regime: in fact, a directly postulated drift-diffusion equation would not capture the exact dynamics in all the possible regimes.

Overall, the numerical results highlight how the competitive or collaborative effects of the cues depend, in a first instance, on the angle between their relative orientations, , between the direction of fiber alignment and the chemotactic gradient. We observed how there is not a simple additivity effect of the two cues and, especially for the case of competitive cues, the determination of the dominant cue strongly depends on their relative strengths, expressed in terms of both concentration and intensity (degree of alignment of the fiber and steepness of the chemotactic gradient). Potentially, the case of competitive cues, combined with the non-local aspects of the model, could lead to interesting further analysis, especially concerning the possible effects of a multi-cue environment on cell adhesion or on collective migration processes.

Eventually, we presented an application of our framework, showing its ability to qualitatively reproduce the experimental results obtained by Sundararaghavan et al. (2013) on competition/collaboration between fibers and chemicals. The presented results suggest the potential applicability of our model to explore further settings, especially thanks to its flexibility in being adapted to the different scenarios. Moreover, its applicability relates to the fact that the presented framework can be used to calculate parameters which quantify directed cell migration (e.g., mean square displacement, persistence time, directional persistence, or mean speed Othmer et al. 1988).

We remark that, even if simulations were performed in a two dimensional setting, the transport models (and their macroscopic limits, as a consequence) are formulated in a general d-dimensional setting. Hence, a possible future development is to perform simulations in the three dimensional case that would be much more realistic for mimicking in-vivo migration of cells in the extracellular matrix. Moreover, our framework is flexible and can be adapted to describe other directional cues that might relate, among others, to haptotactic, durotactic or electrotactic mechanisms. In the same spirit as in Loy and Preziosi (2020), we plan to enrich this framework considering a non-constant sensing radius, as it may vary according to the spatial and directional variability of the external guidance cues. Lastly, this study was restricted to the case in which the cues affect only cell polarization, considering a uniform distribution of the speeds. However, in line with Loy and Preziosi (2020, 2020), this setting may be modified to model a multi-cue environment in which the signals also affect the speed of the cells.