Modelling the collective response of heterogeneous cell populations to stationary gradients and chemical signal relay

The directed motion of cell aggregates toward a chemical source occurs in many relevant biological processes. Understanding the mechanisms that control this complex behavior is of great relevance for our understanding of developmental biological processes and many diseases. In this paper, we consider a self-propelled particle model for the movement of heterogeneous subpopulations of chemically interacting cells towards an imposed stable chemical gradient. Our simulations show explicitly how self-organisation of cell populations (which could lead to engulfment or complete cell segregation) can arise from the heterogeneity of chemotactic responses alone. This new result complements current theoretical and experimental studies that emphasise the role of differential cell–cell adhesion on self-organisation and spatial structure of cellular aggregates. We also investigate how the speed of individual cell aggregations increases with the chemotactic sensitivity of the cells, and decreases with the number of cells inside the aggregates


Introduction
Long-distance collective migration of cells is a biological process important for the establishment and maintenance of multicellular organisms. This orchestrated movement is vital during wound healing, embryonic development and immune responses. Defects during this collective behavior of cells can result in serious health problems, including vascular diseases, tumor formation, and cancer metastasis. These are some of the reasons why the study of the mechanisms by which a population of cells migrate in a coordinated manner is of extreme importance for better understanding many birth defects and diseases [1][2][3].
It is well known that many types of cells respond to chemotactic signals that drive the collective motion of individual cells or cellular aggregates [4][5][6]. Cell chemotaxis has been the subject of intense study during the past years, and thus many details related to the response of individual cells to chemicals are currently known [7]. However, the response of a collective of interacting cells to chemical signals is poorly understood [8]. The complexity of the problem is enhanced by the fact that in many cases these aggregates are com-posed of different subpopulations of cells having different chemotactic responses, and therefore, different motility properties.
In developmental biology, two of the most studied examples of collective chemotactic behavior in different subpopulations of cells are related to the motion of the multicellular microorganism Dictyostelium discoideum, and the motion of neural crest (NC) cells [9][10][11][12]. In regard to the Dictyostelium discoideum cells, their different stages of development are characterized by the cooperative motion of two main types of cells (prestalk and prespore cells) with different chemotactic and mobility properties. Systematic experiments of chemotactic behavior under adjustable and controlled conditions have been performed in order to study the response of groups of Dictyostelium discoideum cells to temporally stable gradients of ′ 3 , ′ 5 -cyclic adenosine monophosphate (cAMP) [13][14][15]. Although many experimental studies removed any chemicals released by the cells themselves, few others considered the interplay between the relay of self-produced cAMP signals and cells response towards stable external cAMP gradients [15]. These studies have been complemented by various theoretical approaches, most of which consider homogeneous groups of cells 2 M Pineda and R Eftimie with similar general properties (e.g. same chemotactic responses and same cell speeds) [14,15]. However, as mentioned above, in many cases collective cell movement is the result of interactions among heterogeneous cell populations with different signal sensing and signal relay capabilities (see Dictyostelium discoideum development that involves prespore and prestalk cells [12]). Similar experiments have been used to analyze the response to chemical gradients of neutrophils [16], cancer cells [17], and bacteria [18].
During embryonic development, different NC cells with different chemotactic responses migrate ventrally through the embryo, guided by migratory pathways [10,11]. Chemotactic interaction between cells is also a very common phenomenon in different physiological contexts. For example, wound healing in the corneal epithelium is determined by the migration into the wound of interacting sub-aggregates of epithelial cells [19], while during tumor progression mutated cells move through and interact with different populations of healthy and immune cells [20].
In this work, we present a general particle model derived to investigate how a mixture of heterogeneous interacting cell populations, which produce, relay and degrade a chemical signal, respond to an external stable gradient of a different chemical signal [21], and form various aggregation patterns. In our model the cells are represented as self-propelled particles (SPP) that behave as soft disks [22,23]. We show how the response of the two populations to the external chemical gradient is determined by the chemotactic sensitivity of the cells, the mutual mechanical interaction between cells, and the overall chemical signal degradation rate [24,25]. Our computational simulations show how self-organisation of cell populations (which could lead to engulfment or complete cell segregation) can arise from the heterogeneity of chemotactic responses alone. These results complement current studies that emphasise the role of differential cell-cell adhesion on self-organisation and spatial structure of cellular aggregates [5,26]. Although this work is inspired by previous studies concerning the response of Dictyostelium discoideum cells to external gradients [15], its novelty lies in the investigation of the interactions between multiple heterogeneous cell populations, and in the investigation of the importance of mechanical versus chemotactic interactions on collective cell migration and aggregation. This paper is organized as follows. In section 2, we introduce the general SPP model implemented in this work. In section 3, we analyze how two groups of mechanical and chemotactic interacting cells, with different chemotactic coefficients and/or average cell speeds, respond to an external stable chemical signal. In particular, we discuss how the relevant parameters of the model affect the response of the cells to different chemical gradient steepness. To this end, we introduce an order parameter that quantifies the global alignment of cells, and we investigate the changes in this order parameter as we vary different control parameters of the model. We also explore how the velocity of single cell aggregations changes with the number of cells and the chemotactic sensitivity. Finally, in section 4 we present our summary and conclusions.

The self-propelled particle model
We describe a group of N cells as a self-propelled particle system, in which each cell behaves as a soft disk of radius = r 0.0075 0 mm. The kinetics of the system occurs in a two-dimensional impermeable × L L x y rectangular domain inside which the state of cell i at time t is characterized by its location vector ( ) → r x y t , , i and by the direction ( ) θ t i of its velocity vector ( ) In agreement with controlled chemotaxis experiments [27,28], the speed ∥ ( )∥ → v t i of each cell is assumed to have a constant value ν i . Then, the updated position of the ith cell is expressed as [7,15,22,23]: , , , and In equation (2) the unit vector → e i gives the direction along which a single cell moves. The first term in equation (3) represents a unitary vector pointing in the direction of increased chemical signal concentration (C T ), multiplied by the chemotaxis sensitivity χ i . Normally, it is assumed that χ i decreases nonlinearly with C T [29][30][31]. However, according to the experimental results that motivated this work [15], we will assume that this term is independent of the level of the chemical signal. Note that arg gives the angle of the resulting vector.
In our model the total chemical signal C T is originated from two sources: the single cells that secrete a local chemical signal C I into the extracellular space, and an externally imposed stable gradient of a different chemical C E which is changing only along the y direction. We consider that After assuming that secretion, degradation, and diffusion of C I are much faster than all other processes of the model, the chemical gradient generated by the cells can be written as where K 1 is a modified Bessel function of the second kind, / µ = D I I , and p quantifies the response of the cells to the chemical gradient. The sum is over the simulated cells other than cell i. The vector → r ij is a unitary vector directed from j to i. The parameters D I and µ I are the constant chemical diffusivity and chemical signal degradation rate, respectively [32]. Note that the parameter N is the total number of cells and the local chemical signal is secreted at every position r i of cells.
We assume that the external chemical signal only changes along the y directions and its dynamics is governed by with boundary conditions along the y direction: As in [15], we choose fixed-concentration boundary conditions to describe the case where the external chemical is applied only at the upper boundary of the domain. Assuming also steady state conditions (since the dynamics of the chemical is much faster than all other processes in the model), one obtains with the corresponding gradient given by( where ŷ y is a unitary vector point in the y direction. Note that the gradient of the external chemical described in equation (7) depends on the concentration C m at the upper boundary of the domain. In supplemental material we plot equation (6) (external chemical profile) as a function of L y (figures S1(A)) and µ E (S1(B)). Figure S1(A) shows that the chemical profile substantially increases near the upper edge of the domain where ( ) = = C y L C E y m . However, figure S1(B) reveals that the steepness of such a profile decreases with the chemical degradation rate µ E . In particular, it shows that for small enough values of the degradation rate the external chemical profile may extend along the whole domain.
The second term in equation (3) is responsible for the repulsive soft mutual interaction between cells [33,34]. To describe it, let us consider a cell i. A second cell j exerts a repulsive force ( ) , , ij i j is the distance between cells. The repulsive force decreases linearly with distance from the center of the cell, and it is zero for ( ) > r t r 2 ij o (where we define r 2 o to be the interaction radius between two cells). The parameter m ij that appears in equation (3) controls the intensity of the repulsive cell-cell interactions.
The last term in equation (3) is described by , a random uniformly oriented unitary vector, multiplied by parameter α i , which quantifies the noise intensity. The variable ( ) ζ t i is a random number drawn between 0 and π 2 . Finally, in this work we assume for simplicity that the chemical diffusivity and degradation rate are the same for both the local and external chemicals: However, the response of the cells to the local chemical signal is distinguished from the response to the external chemical gradient by the parameter p. We also assume that = m m ij c and α α = ∀i j , i . With these assumptions in mind, we present the results of our study in the following section.

Results
As mentioned above, in this work we investigate how heterogeneous interacting cell subpopulations with distinct motility properties react towards externally imposed chemical signals. To achieve this, we divide the group of N cells into subpopulations A and B [24,25]. For simplicity, let us also assume that both subpopulation have the same number of cells ( ), and that cells inside each subpopulation exhibit the same mechanical and chemotactic properties. However, given that the subpopulations may differ in cell speed and/or chemotactic sensitivity, we consider each cell i of type A or B (where cells A have lower speeds and/or chemotactic sensitivities compared to cells B). Then, the system is iterated in time with a time step ∆ = × − t 2 10 2 min. To be in line with the experiments of Guven et al [15] in which a uniform cell injection flux is considered (with the orientation of each newly introduced cell assumed to be initially in the y direction), we simulated equation (1) with all cells starting at = y 0, uniformly distributed along x axis and oriented only in the y direction. However, after this initial step, the cells can change their orientation (in response to the chemical signal, the social forces, and randomly). Then, we proceed to characterize the collective motion of cells by the order parameter M defined as in [35][36][37]: In the following section we start exploring the collective alignment of the model when the two groups differ only in their chemotactic sensitivity. Then, we consider the case in which the groups differ in chemotactic sensitivity and cell speed. In the last subsection, we investigate how the speed of a single cell aggregation changes with relevant control parameters.

Two groups of cells with different chemotactic sensitivity
In this section we consider the case in which the two groups of cells exhibit different chemotactic sensitivities: each cell i belonging to the group A reacts to the chemical signals with a chemoattractive sensitivity of χ χ = = 0.75   but also a complete spatial sorting of the two subpopulations of cells emerge (with the majority of the less chemotactic aggregates lagging behind the strong chemotactic ones). In panel (C) we see that slightly larger values of L y lead to undirected and engulfed aggregates (the stronger chemotactic cells are engulfed by the less chemotactic ones). Thus, increasing the size of the domain along the y direction promotes the formation of undirected cell aggregations each of them containing both populations. Figure 2 shows that the extension of the domain along the x direction (parameter L x in our model) increases the number of cell aggregations. It also suggests that small and strong chemotactic cell aggregations move faster (we will come back to this point later). Note that in this case, simulation conditions and parameters are as in figure 1. From figures 1 and 2 one can conclude that, under our simulation conditions (i.e μ and D fixed), the larger the size of the domain the larger the number of undirected cell aggregations formed. Now we proceed to study how the alignment of the cells is affected by the control parameters p, μ, and C m . The spatial cell configurations presented above emerge from the interplay between the cell-cell  We would like to emphasize that the merging process between aggregates is more frequent when the influence of the external chemical gradient is small. Under this condition, the aggregates have enough time to interact and merge. When the influence of the external gradient is large, the cells tend to move aligned towards the y directions and the aggregates do not have time to merge before reaching the upper wall of the domain. We note that the aggregates also merge when the size of the domain and/or the cell-cell attractive chemical interactions are large enough.  aggregate can emerge and remains inside the domain for very long time-shown here is = t 400 min. This single aggregate is the result of the collision between two separate aggregates due to the mutual chemical attraction (see the corresponding two aggregates in figure 4(C)).
One can associate the black and dark blue regions of figures 3(A') and (B') with the formation of slow moving aggregates, the yellow regions with a more or less independent migration of the cells towards the external gradient, and the transition regions with the formation of aggregates that most of the time move as a whole towards the external gradient. Aggregate formation was assessed by visual inspection.
In the following section we explore the case in which the two sub-populations also differ in their average cell speeds.

Two groups of cells differing in cell speed and chemotactic sensitivity
Taking into account the common assumption that strongly chemotactic cells move faster [6], we assume that each cell i belonging to the group A reacts to the chemical signal with a chemoattractive sensitivity Here we show the case in which the two groups differ in velocity and chemical sensitivity (corresponding to figure 6(C)). In this case the black arrows represent slow and less chemotactic cells, while the red arrows represent fast and highly chemotactic ones. . First, let us mention that the impact of the size of the domain and the parameters p, μ, and C m on this new version of the model is similar to the one reported in previous section. Namely, the global alignment of cells decreases with the size of the domain along the y direction, the chemo-attraction between cells, and the chemical degradation rate, and increases with the steepness of the external chemical signal. Therefore, in the following we report only the major differences between the two cases (with similar or different cell speeds).
In figure 6 we plot separately the temporal evolution of the order parameter M for each subgroup normalized by the total number of cells (i.e. Our simulations show that, in the long run, the two group of cells align almost totally. This is in agreement with figures 3(B) and (B') of the previous section. The numerical simulations also show that there is a time interval during which the slow and less chemotactic group of cells exhibits a better alignment than the fast and stronger chemotactic group of cells. Interestingly, the period of time during which this phenomenon occurs increases with the number of cells in the system. However, as figures 6(D)-(F) reveal, this phenomenon does not exist or it is not significant when the two groups differ only in their chemotactic sensitivity. In this case, after a short time the two groups exhibit almost the same alignment. Figures 6(G)-(I) show that this interesting phenomenon occurs even when the two groups have the same chemotactic sensitivity, but differ in their average cell speed. In both cases, the duration of the migrating phenomenon increases with the number of cells in the system. In figure 7, we show the spatial cell distribution of = N 500 cells, for the case in which the two groups differ in their average cell speed and their chemical sensitivities (corresponding to figure 6(C)). It shows that at = t 70 min the fast and stronger chemotactic red (gray on black/white prints) cells collect in the middle of the aggregates, with their velocity vector pointing in random directions. In con-trast, the slow and less chemotactic black cells move aligned towards the aggregate formed by the inner cells or/and the external chemical gradient (see videos S3 in the supplementary material, for a better understanding of the corresponding cell distribution patterns). Finally, we would like to emphasize that for all these plots the time runs from 0.02 min to 100 min. This is the reason why the initial state in which = M 1 does not appear in these simulations.
We also analyze the impact of the noise parameter α on the cell migrating phenomenon presented above. Figure 8 shows the temporal evolution of M for three different noise intensities, when = N 500. Note that the time interval during which this phenomenon occurs increases with the noise intensity. This is due to the fact that the noise delays the time for cells to reach a completely aligned state, as induced by the external chemical gradient. These results also reveal that, as expected, the difference in the levels of alignment between the two groups decreases with the noise intensity.
In the following section we explore the changes in the speed of cell aggregations as a function of some model parameters.

Speed of cell aggregations as a function of control parameters
In this section we explore the speed of single aggregates as a function of the number of cells (N) and some other parameters of the system (e.g. p, μ, L y ). For simplicity, we assume that all cells have the same average speed (ν , for all i). We pay attention to the velocity component along the y direction of a single aggregate (V y ). To calculate this velocity, we identify the center of the aggregate and measure the displacement of this central point during the last quarter of the simulation (150 min < < t 200 min). We assume that initially, cells are homogeneously distributed in a small sub-domain, with all of them pointing in the y direction. We denote the number of cells in the aggregate by N (since for these simulations all cells in the population form an aggregate). In the supplementary material we show examples of some aggregates used to calculate the speed (e.g. see video S4 for an illustration of how these aggregates move over time and space).
In figure 9 we plot V y as a function of the cluster size N, for several values of p, μ, and L y . The figure shows that the velocity along the y direction decreases with the number of cells (since the attraction towards the center of the aggregate due to local chemical signals increases with the increase in cell numbers). Panels (A), (B), and (C) also show that such velocity decreases with the response of the cell to the chemical signal secreted by themselves (parameter p), the chemical degradation rate μ, and the extension of the domain along the y direction, L y . Moreover, figure 10 shows that V y increases with the chemotactic sensitivity of cells. Also in this case, the speed decreases with p and μ.
These results explain why, in the case of two different subpopulations, the smaller and stronger chemotactic aggregates move faster (see for instance figure 2).

Summary and discussion
In this work, we generalized a one-population selfpropelled particle system introduced in [15] to further describe systematically the dynamics of populations of distinct heterogeneous cells, and then used it to assess the chemotactic response of the cells to imposed stationary chemical gradients. First, we considered two cell populations that differ in their chemotactic sensitivity. Then, we analyzed the case in which the two populations also differ in their average cell speed. We found that the response of the two populations to the external chemical gradient was strongly affected by the response to the chemical signal secreted by the  We also found that the collective alignment of cells decreased with: the size of the domain along the y direction, the response of the cells to the chemical signal generated by themselves, and the chemical signal degradation rate (see figures 1 and 3). Our simulations also showed that small steepness of the external chemical gradient favours the formation of cell aggregates and the undirected motion of cells, while large gradient steepness induced complete cell alignment (see figure 3). Moreover, the simulations commonly showed the formation of cell aggregates in which the fast moving cells were engulfed or partially engulfed by the slow moving ones. In addition, a number of completely segregated cellular aggregations were also observed (see, for instance, figures 1(B), 2, and 4, and the videos in the supplementary material). We also found that the merging process between cell aggregates was more pronounced when the influence of the external gradient was not too strong to counteract the chemical aggregate-aggregate attraction (see figure 5).
Although not discussed in detail, we also found that the phenomena mentioned above occurred in the situation where the two subpopulations of cells simultaneously had different chemotactic sensitivities and cell speeds (see figure 6 for an example of complete alignment that both versions of the model exhibited in the long run). However, in this work we only payed attention to the main difference between the two aforementioned cases. We found that when the two groups differ in cell speeds and chemotactic sensitivities, there was an initial-to-intermediate period of time during which the slow moving cells exhibited a better alignment than the fast moving ones (see figure 6). This behaviour, which was not observed when the subpopulations differ only in chemotactic sensitivity, was mainly due to the differences in cell speed between the two subgroups of cells. We also found that the intrinsic noisy behaviors of cells lead to an increment in the time interval during which this phenomenon was observed. However, the noise also decreased the difference in alignment between the two subpopulations (see figure 8).
Finally, we analyzed the speed of a cell aggregate as a function of the main control parameters of the model. For simplicity, we assumed that all cells possessed the same average speed and chemotactic sensitivity. We generated a single aggregate and measured the velocity of it while moving along the y direction towards the external chemical gradient. We found that the velocity decreased with: the number of cells, the response of the cells to the local chemical signal generated by themselves, the chemical degradation rate, and the extension of the domain along the y direction (see figure 9). Our simulations also showed that the speed of aggregates increased with the chemotactic sensitivity of cells χ (see figure 10). These results allow us to conclude that, in our model, smaller and highly chemotactic cell aggregates move faster compared to larger or weaker chemotactic aggregates.
The model has a number of experimentally relevant parameters (e.g. N, p, μ, ...), which allowed us study numerically, in a simple but general way, the collective migration of heterogeneous cell populations. This study was motivated by the experimental results reported by Guven et al [15], which investigated the response of a single homogeneous population of cells to external chemical linear gradients and signal relay. Thus, in principle, we expect that some of the phenomena reported in our study also occur in experiments performed with subpopulation of cells having different chemotactic properties and/or cell speeds (see, for example, the chemotactic cell sorting reported in [1] for a heterogeneous population of skin fibroblasts and malignant fibrosarcoma cells). We also need to emphasize that our model does not take into account the exact experimental conditions implemented by Guven et al [15]. In that exper imental study the authors injected a uniform flux of cells, in which any newly introduced cell was oriented more or less in the y-direction. Thus, the number of cells in their experimental setup increased in time, which is in contrast to our theoretical model where we assumed a finite number of cells at the initial time. Furthermore, in Guven et al [15], once the cells reached the upper wall of the domain, they were removed from the experiment. That experiment also considered that the external chemical gradient is a linear one, whereas in our theoretical model the gradient normally exhibits an exponential decay. All these differences could be included in our minimalistic model. We will explore these aspects in our future research.
Although the model investigated in this study predicted interesting and potentially verifiable collective phenomena in heterogeneous cell populations (as is the segregation of stronger/weaker chemotactic cells in [1]), it could be further extended in several different ways. Two such approaches, which focus on the incorporation of saturated cell chemotactic responses to large concentrations of chemicals (as observed experimentally [13]), and on the incorporation of cell-cell adhesive forces (which have been shown to affect cell aggregation and sorting [5,38]), are also the topic of on-going work. Instead of considering two single subpopulations of cells, it could be also interesting to assume that each cell has a different response to the chemical signals. These type of cell-to-cell variations were recently included in a model of collective chemotaxis, and it was found that chemotaxis is limited by cell-to-cell variation in signaling [39]. While our study focuses on a simple model, our numerical results could help advance the understanding of the complex and orchestrated movement of heterogeneous cell populations toward external chemical signals.