Enhancement of Chemotactic Cell Aggregation by Haptotactic Cell-To-Cell Interaction

The crawling of biological cell is a complex phenomenon involving various biochemical and mechanical processes. Some of these processes are intrinsic to individual cells, while others pertain to cell-to-cell interactions and to their responses to extrinsically imposed cues. Here, we report an interesting aggregation dynamics of mathematical model cells, when they perform chemotaxis in response to an externally imposed global chemical gradient while they influence each other through a haptotaxis-mediated social interaction, which confers intriguing trail patterns. In the absence of the cell-to-cell interaction, the equilibrium population density profile fits well to that of a simple Keller-Segal population dynamic model, in which a chemotactic current density J→chemo∼∇p competes with a normal diffusive current density J→diff∼∇ρ, where p and ρ refer to the concentration of chemoattractant and population density, respectively. We find that the cell-to-cell interaction confers a far more compact aggregation resulting in a much higher peak equilibrium cell density. The mathematical model system is applicable to many biological systems such as swarming microglia and neutrophils or accumulating ants towards a localized food source.


Introduction
Understanding the mechanisms behind cell population dynamics is essential to a wide range of biological processes including development, wound healing, tumor expansion, and immune responses [1][2][3][4][5][6][7]. However, it is a very challenging task, since the relevant systems generally involve many different interacting constituents (cells) which are inherently nonlinear and the type of cell-to-cell interactions vary from one case to another. For example, neighboring cells can interact with each other by mechanical forces in a close-packed layer of cells [8,9]. The local forces can then generate long-range spatio-temporal correlations. Some well-known examples include the ratchet-like tissue movement during dosal closure in the developing Drosophila embryo [10] and the waves and swirls in the in vitro systems of an expanding epithelial cell sheet [2,11,12]. Cells can also be coupled through diffusing chemical agents and matching receptors. One of the most well-studied examples of this type of cell is the traveling-wave chemotaxis of dictyostelium discodium (or dicty) amoebae [13][14][15]. Briefly, the starvation triggers amoebae to produce and excrete 3',5'-cyclic adenosine monophosphate (cAMP) that diffuses to the neighboring cells which have cAMP receptors. Not only can the cells amplify the level of cAMP, they can also dissociate cAMP to cGMP in a temporally coordinated manner. Consequently, the cAMP-mediated coupling can bring about large-scale cAMP waves. Therefore, for the case of dicty amoebae, the diffusive coupling and the cell-intrinsic nonlinear kinetics and adaptation are responsible for the collective phenomenon. The amoebae cells also actively move (i. e. chemotaxis) towards the higher concentration of cAMP while experiencing the positive slopes of cAMP waves.
In some cases, the chemical agents released by crawling cells do not diffuse but stay behind the cells, and become encapsuled into small vesicles (typical size, 40*100 nm in diameter) called as exsosomes. Alternatively, they become bound to the two-dimensional substrate (or three-dimensional matrix) on which the cells are placed [16][17][18]. In our recent work, we showed that in cultured microglia, the immune cells of the brain [19], a number of chemical markers (e.g. α5-integrin) are left behind the crawling traces of the cells and that these non-diffusing (but degrading) markers serve as intercellular signaling agents guiding their motility [20], enabling the cells to form intricate evolving network of trails [21]. In short, microglia exhibit "haptotaxis" by following the trails they generate. Generally speaking, "haptotaxis" refers to the directed movement of cells controlled by the relative strengths of peripheral adhesions into some substrates.
Microglia also exhibit "chemotaxis" in response to the concentration gradient of adenosine monophosphate (ATP) [22,23]. "Chemotaxis" is the directed movement of bacteria, eukaryotic cells, or multi-cellular organisms toward concentrations of environmental chemoattractants. Therefore, microglia movements can be driven by several different agents. Multi-modal communications are indeed very common for biological cells. Neutrophils rely on integrins and lipid leukotriene B4 (LTB4) for swarming [24]; they also release exsosomes behind their moving trails as for intercellular signaling [16]. Numerous chemokines are known to be involved in chemotaxing tumor cells [7] and tumor cell migration is driven by haptotaxis as well as chemotaxis [25]. Indeed, for many realistic biological settings, cell-to-cell interactions are mediated by several different modalities with many different intercellular signaling molecules working together. However, the potential interplay between these different modalities is largely unexplored. This paper, in particular, investigates the role of a haptotaxis-mediated cell-to-cell interaction in the aggregation dynamics of the chemotaxing of the cells in response to a globally imposed chemoattractant gradient. The paper is based on the computer-simulations of a population of model cells having some key features of crawling microglia in in vitro culture. More specifically, this work is motivated by our recent experimental observations concerning microglia trail network formation and the chemotactic behavior of microglia. While the haptotaxismediated trail network formation and the chemotactic responses of microglia to ATP have been studied separately in the past, their interplay has not yet been addressed. Thus, this paper may be considered as an exposition of the ways in which such an interplay can be interesting and useful. Most importantly, we find that the haptotaxis-mediated trail network can greatly enhance the aggregation of chemotaxing cells.

Results and Discussion
Properties of the rule-based mathematical model cell mimicking a microglial cell surface (see Fig 1a). The movement of the loop is determined by two scalar fields: the activator S þ ðr; tÞ and the inhibitor S − (t). S þ ðr; tÞ depends on the positionr (the position vector along the cell perimeter with respect to the centroid of the cell) as well as time, whereas S − (t) is a global variable that only depends on the time. At every iteration time step,r (along the perimeter) can either advance, retreat or stay based on the following set of rules. Retraction occurs when S þ ðr; tÞ S À ðtÞ and the rate of retraction is governed by the following equation: where r min is the constant minimum radius and R − is the retraction rate constant. The function max(x, y) selects the larger value of x and y. Protrusion occurs when S þ ðr; tÞ > S À ðtÞ at a rate governed by the following equation: where R + is the average protrusion rate and G(R + ) is a random number generated from a Gaussian distribution with the mean R + and variance R + . The activator S þ ðr; tÞ evolves according to the following equation: @S þ ðr; tÞ=@t ¼ K diff r 2 S þ ðr; tÞ À K decay S þ ðr; tÞ þmax½Gðf ðS þ ðr; tÞ À S À ðtÞ; g; l þ P baseline ÞN burst Þ; 0: The first term accounts for the diffusion of S + along the cell perimeter and the second term renders a self-decomposition of S + . The last term is a stochastic positive feedback loop accounting for both the local stimulation and the existence of a random signal. The function f(x, γ, λ) = 0 for x < λ and (x − λ) for x ! λ, where λ is a threshold value for the feedback. P baseline refers to the rate of random bursts caused by the internal baseline activities. The function G again represents a random number generated from a Gaussian distribution.
The retraction signal is governed by the global inhibition rule S À ðtÞ ¼ C À A R S þ ðr; tÞdr, where C − is the inhibition constant, A is the total area of the cell, and the integration is a line integral over the entire cell border, which is composed of 360 pixels (i.e. 1 pixel for 1 degree with respect to the centroid). Each pixel corresponds to 0.286 μm and one iteration time step is 1 s. At each iterative time step, the formation of focal adhesions and their detachments are assigned stochastically to the points along the cell perimeter with a probability P þ fa and P À fa , respectively. Retraction is inhibited when a perimeter point hits a focal adhesion. Again, all the rationales for the kinetics and mechanics are fully described in the original paper by Satulovsky et al. [26].
It is known that the rule-based phenomenological model can successfully recapitulate many different types of cell morphologies and motilities, mimicking cells such as keratocytes, neurons, and amoebae. In the current study, the following set of parameter values is chosen to match some of the important characteristics of a typical microglia: R + = 0.103 μm/s, R − = 0.0281 1/s, K diff = 11.9 μm 2 /s, N burst = 13, P baseline = 0.181 1/(sÁ μm), λ = 3.22 1/μm, γ = 29.1 1/s, C − = 1.93 × 10 −5 1/μm 3 , P þ fa ¼ 0:0003, P À fa ¼ 0:0058, and r min = 2.857 μm. Only the parameter K decay is varied as a control parameter. Fig 1a shows two exemplary passages traced out by two freely crawling model cells having different values of K decay . As clearly illustrated in Fig 1b, for the time domain ranging from hundreds to a few thousand seconds (yellow shaded area in Fig 1b) the cells are superdiffusive with the exponent γ of the mean squared displacement <R 2 (t)> * t γ with the range of 1 < γ < 2 (see Fig 1c, green dots). Overall, γ decreases as K decay increases, and this tendency is reflected well in the two exemplary traces shown in Fig 1a. The directional persistence of a given passage can also be estimated from < cosθ(t) > vs. t, where θ(t) is the angle between two tangent vectors obtained at two different locations along the passage that are separated by time t (see Fig 1d): < cosθ(t) > is fitted to a function Ae Àt=t 1 þ ð1 À AÞe Àt=t 2 ; and then τ 2 is then considered as the persistent time of crawling, while τ 1 is an estimate of the small time interval (typically, 3*7 minutes) between two successive small-angle zigzag turns (see Ref. [27] for details). Similar to the superdiffusiveness exponent γ, overall the persistent time τ 2 is also a decreasing function of K decay , as shown in Fig 1c. Chemotaxing model cell population under a globally imposed concentration gradient of an attractant Motivated by the earlier experimental studies showing chemotactic behaviors of microglia towards a high concentration of ATP [22,23], we modified the single cell model to be responsive to an extrinsically imposed chemotactic drive; specifically, we added a new term pðr 0 Þ to Eq (3), wherer 0 represents a two-dimensional position vector with respect to a laboratory frame and pðr 0 Þ is a concentration map of the chemoattractant. The modified equation for S + reads, @S þ ðr; tÞ=@t ¼ K diff r 2 S þ ðr; tÞ À K decay S þ ðr; tÞ þmax½Gðf ðS þ ðr; tÞ À S À ðtÞ; g; l þ P baseline ÞN burst Þ; 0 þ pðr 0 Þ: As the level of p increases, the production of activator S + also increases. Consequently, the cell boundary facing upwards of the gradient (i. e. higher value of p) is more likely to have a higher value of S + than that of the cell boundary facing downwards of the gradient. On the other hand, S − is a global variable depending on the spatial average of S + . Therefore, it will not be greatly affected by the p gradient. Subsequently, the cells are more likely to move up the gradient on the average.
In Fig 2a, we placed 2,000 non-interacting model cells uniformly but randomly, and imposed a conically shaped chemoattractant pðr 0 Þ concentration image [the opening angle of the cone, θ = 2 arctan (1.25r 0 )] to the system to investigate its subsequent chemotactic response. For the given p gradient, although their long-term motilities are biased towards the peak of p on average, the individual model cells still trace out quite random passages locally due to their built-in stochastic elements: An exemplary trace is given in Fig 2b. In the long run, the population cell density σ equilibrates with a peak in the middle and it fits fairly well to an exponential function as shown in Fig 2c (see S1 Movie). To better understand the origin of this equilibrium profile, we measured the chemotactic (surface) current density K chemo as a function of σ and the diffusive current density K diff as a function of @s @r 0 for several different values of K decay as shown in Fig 2d and 2e, respectively. Here, K chemo (or K diff ) is measured to be the number of cells, moving up (or down) the p-gradient, crossing a unit arc-length along a circle about the peak per unit time.
The measured K chemo is a linear function of σ. This seems quite natural, because the cells are essentially identical; thus, they have the same crawling tendency towards the peak of p and for the given setting there is no cell-to-cell interaction. Therefore, K chemo should be a linear function of σ. Indeed, the linear form K chemo = ασ (with α being the chemotaxis strength) has been implemented widely as the simplest form of chemotactic current density [28]. The chemotatic strength α, which can also be viewed as a net chemotatic velocity, is a slowly decaying function of K decay (see the inset of Fig 2d). Considering that the model cells are built around nonlinear kinetics, rule-based nonlinear mechanics and several stochastic features, the nearly linear relationship cannot be explained easily. Surely, α will depend on many different factors such as instantaneous crawling speed, directional persistence in crawling, angular diffusivity, and memory.
On the other hand, K diff is measured to be a linear function of the density gradient @s @r 0 as shown in Fig 2d. The slope of K diff vs @s @r 0 is nothing but effective diffusion coefficient D; it is evaluated for several different values of K decay and is shown in Fig 2f (blue line). In other words, although the individual model cells are superdiffusive and have a directional persistence for the  Fig 2e (Fig 1b). The red dashed line in (c) represents an exponential fit of the data. The variations in α and D brought about by a different initial condition are less than 0.5% (n = 5). Fig 1b), for the long time domain (t ≳ 10,000 s, shaded in violet in Fig 1b), which is relevant for establishing the density equilibrium state, they can be viewed simply as a Brownian particle following the normal Fick's diffusion law. In fact, similar values of D are also obtained by fitting the long time domain portions (10,000 < t < 20,000 s) of the MSD curves of Fig 1b to a straight line with slope 1 as shown in Fig 2f (red  line). Within the estimated error bars associated with the linear least-square fitting, the two different estimations agree fairly closely. D is measured to be a decreasing function of K decay or an increasing function of the directional crawling persistence. Within the long time regime, in which the cells can be viewed simply as a random-walker, a decreased directional persistence in the crawling cell is effectively similar to having a smaller step size for a random walker. Besides, we found that the angular diffusivity increases as K decay increases (not shown, see Ref. [27]). Consequently, parameter K decay strongly influences D and Fig 1a also clearly illustrates this point.

Trail formation by haptotaxing model cell population
We previously showed that crawling microglia can form intricate networks of trails (see S1 Fig) and the mathematical model cells can faithfully recapitulate many of the key features of microglial trail networks if they are coupled by a haptotaxis-mediated coupling [21]. Here, we reproduced the phenomenon in a disk geometry (see Fig 3). For the haptotaxis-mediated coupling, we added a chemoattractant to Eq (3). qðr 0 ; tÞ is a scalar field representing the concentration of a number of proteins produced and deposited atr 0 by crawling cells. t i represents the specific time of ith visit to the positionr 0 previously made by one of the cells. After being deposited, the attractant selfdecomposes with a time constant τ decay . From the viewpoint of S + kinetics, the role of q is essentially identical to that of the chemoattractant p: When a model cell encounters a q-deposited trail, it will most likely be guided along the trail since q is a S + enhancing chemical agent similar to p. The scalar field q is, however, a dynamic variable that evolves spatio-temporally, while p represents the extrinsically imposed fixed landscape of a chemoattractant. Again, from an initially uniform, randomly-distributed state (Fig 3a), the model cells are allowed to move and release q behind their passages. As they crawl around and recognize their neighboring q-landscape, their motile behaviors become biased to follow the nearby "foot prints" created by other cells. Consequently, an evolving network of trails is formed, as illustrated in Fig 3b. For the chosen set of parameter values, small trails initially emerge, forming a dense fuzzy network (Fig 3b, t = 3 hr), which eventually becomes a more refined network of trails (Fig 3b, t = 24 hr, also see S2 Movie). The network evolves constantly as new branches form, while some less visited portions disappear (compare the two pairs of boxed areas in Fig  3b). From a uniformly distributed state, the establishment of a mature network typically takes about 24 hrs, as shown in Fig 3c, in which the trails are defined to be the area with q > 0.22. Subsequently, the trails are "skeletonized" with an open source program ImageJ and each crossing junction is identified as a node.   values of K decay (see Fig 4c), the skewness changes significantly as expected (see Fig 4d). The mean instantaneous crawling speed along the trails is almost a linearly decreasing function of K decay , ranging approximately from 5.7 to 6.7 μm/min: As the trail becomes straighter, the speed at which the cells move increases. However, the speed at which the model cells crawl is about 30% faster in the absence of any q-trails (see S2 Fig).
Effective aggregation assisted by cell-to-cell interaction resulting trail network Finally, we question the consequence of implementing the haptotactic cell-to-cell interaction for the cells that are already chemotaxing. The cells experience the same static p-gradient as shown in Fig 3a (see Fig 5a, t = 0), but they now also have the q-mediated interaction forming trails. Due to the chemotactic force of the p-gradient, initially the self-organized trails exhibit some tendency to orient themselves towards the center (see Fig 5a, t = 6 hr, also see S3 Movie). On the other hand, the presence of the trail network can break the circular symmetry of the system and render the centroid of the "high density core area" not in the middle of, but offcenter to, the peak of p. The core-centric chemotactic force competes with the non-centric haptotactic guide imposed by the trails. For the chosen set of parameter values and the plandscape, we find that the high density core centroid is typically 32%*35% off-centered while its azimuthal location changes with a different initial state of the cell population (see S3a Fig). Importantly, however, the position of the (time-averaged) maximum cell density is still located at the peak of p (see S3b Fig), which is also the global maximum of p + q (see S3c Fig).
Fig 5b compares two equilibrium density profiles as a function of the radial distance from the peak of the population density, which is also the center of the disk domain: one corresponds to the case having the q-mediated haptotactic interaction (red line) and the other corresponds to the case without it (blue line). Clearly, the cells are aggregated more around the peak (r 0 = 0) for the case with the q-mediated social interaction. In fact, in equilibrium, the number of cells in the core area (bounded by the dashed circle in Fig 5a 24hr) is almost two times that of the case with no q-mediated interaction (see Fig 5c). In other words, the existence of trails markedly improves the cell gathering in the core area. This can be attributed to the structure of the trail network (see S3b Fig for example) that encircles the peak core area. The structure effectively reduces the diffusive flux down the hill, guiding the cells along the azimuthal direction. On the other hand, the chemotactic uphill flux is not so much affected by the existence of the trails as the cells channels through the trails towards the peak. In equilibrium, the majority of cells are being trapped within the core area, where the level of p + q is very high (see S3c Fig).

Conclusions
We conducted numerical simulations of the population dynamics of mathematical model cells that have a built-in directional persistence in their crawling for three different cases: 1) when they perform chemotaxis following an externally imposed gradient of a chemical attractant, 2) when there is a haptotaxis-mediated cell-to-cell social interaction, and finally 3) when both are present simultaneously. Originally, the model system was developed to understand and mimic the properties of microglia that exhibit both chemotaxis and haptotaxis. However, while it seems that all cell movements are generally driven by several different factors and mechanisms working simultaneously, very little has been discussed so far about their collaborative interplays. This work may be considered as an exposition of the way in which such an interplay can be synergistically useful or important in biology. Most significantly, we demonstrated that a haptotaxis-mediated cell-to-cell interaction could significantly enhance the chemotactic aggregation of cells towards the area of maximum concentration of a globally imposed attractant. Several other results were also notable. For the case of 1) (chemotaxis only), the steady-state cell density profile fits rather well to that of the simple Keller-Segel chemotaxis model [28], in which a chemotactic flux J chemo = αρ competes with a diffusive flux J diff = Drρ to result in a density profile of an exponential function. For the slow time scale which is relevant for reaching an equilibrium density profile, the individual model cells can be viewed simply as a normal Brownian particle. Considering that no cell-to-cell interaction occurs, the linear dependence of J chemo on ρ seems natural. It is also worth noting that in the case of chemotaxis with a haptotactic social interaction, the aggregation core (i.e. the centroid of the population) is in general not centered on the peak of the chemoattractant p profile and the overall azimuthal symmetry is broken. Therefore, a symmetry-breaking instability seems to be responsible for producing a localized density structure on the trail network by a self-enhancing feedback mechanism. This needs to be further explored.