Modeling chemotaxis of adhesive cells: stochastic lattice approach and continuum description

The effect of chemotaxis on migration of adhesive and proliferative cells on a substrate is analyzed by employing two approaches: by solving a stochastic discrete lattice model for cell dynamics and by deriving and solving a continuum macroscopic equation for cell density. The phenomenon of front propagation is investigated in the framework of the two approaches both for positive and negative chemotaxis. A good agreement between the results of the lattice model and of the continuum model is observed both for front velocities and front profiles. The theoretical model is also able to match recent experimental observations on glioma cell migration.


Introduction
Biological multicellular systems present an exciting example of active matter, which is intrinsically out of equilibrium. Living cells or bacteria can migrate individually nevertheless, producing the whole spectrum of emergent collective behavior ranging from snow-flake like bacterial patterns [1,2] to aggregation of myxobacteria [3,4] and clustering of brain tumor cells [5]. Glioma (brain tumor) cells are highly motile. If placed on a substrate, a cell can migrate its own diameter (order of 20 µm) in 3-10 min [6]. This ability of cells to migrate and invade the extracellular matrix moving away from the main tumor mass makes brain tumors very difficult to treat [7,8]. 1 Author to whom any correspondence should be addressed.
In addition to random motion, a cell can perform directed motion in response to a gradient of some chemicals; this phenomenon is called chemotaxis [9]. Chemotaxis can be positive if a cell moves toward a higher concentration of chemoattractants; during negative chemotaxis cells move toward the lower concentrations of chemorepellents. Bacterial chemotaxis is an active area of research [10]. The ameba Dictyostelium discoideum employs chemotaxis for intercellular interactions. Dicty secretes cyclic AMP attracting other bacteria since they migrate toward a higher concentration of cAMP [11].
Glioma cells also exhibit chemotactic behavior, moving toward higher concentrations of various growth factors or glucose. It was shown that gliomas are chemotactically attracted to epidermal growth factor (EGF) without any effect of EGF on cell motility (random cell motion) [12]. Another growth factor that leads to directed migration of glioma cells is T G Fβ 1 [13]. In a series of works Del Maestro and co-workers [14,15] experimentally analyzed migration of glioma cells away from a tumor spheroid introduced into gel. The authors found evidence that glioma cells secrete a chemorepellent factor; in particular, when two tumor spheroids were introduced, cell invasion into the region between the spheroids was substantially smaller, possibly due to higher concentration of chemorepellent there [15]. On the modeling side, one can employ stochastic lattice models for cell migration, describing negative chemotaxis by introducing effective repulsion between cells [16]. Another approach is to formulate a system of reaction-diffusion equations for cell density and for concentration of chemoattractant (or chemorepellent) factors [17,18]. Finally, one can employ a hybrid approach, treating cells as discrete objects migrating on a lattice, and formulating a continuum equation for chemicals [19].
In this work we analyze the phenomenon of cell migration (in terms of front propagation) by deriving a continuum equation for cell density from the underlying microscopic lattice model taking into account both cell-cell adhesion and chemotaxis. The theoretical results obtained by solving the system of reaction-diffusion equations agree very well with the numerical results of the stochastic hybrid model. In a recent work we investigated glioma cell migration on a substrate, both experimentally and theoretically [20]. The theoretical model that incorporated cell motility, cell proliferation and cell-cell adhesion was not able to produce results that matched the experimental observations. Rather, cell density would drop off too quickly compared to experiment. Here we show that when cell chemotaxis is taken into account, the theoretical results agree with the experimental data.

The discrete lattice model
To describe cell migration on a substrate, we employ a stochastic lattice model that incorporates cell diffusion (random walk), cell proliferation, cell-cell adhesion and chemotaxis. Cells are placed on a two-dimensional square grid on which they can move up, down, left or right. Every lattice site can be empty (n = 0) or occupied by a single cell (n = 1), incorporating the exclusion principle. A cell performs random walk on a lattice with a characteristic diffusion time t diff of the order of 5 min, however only jumps to empty neighboring sites are allowed. We define up/down as the y-axis and left/right as x. Cells are initially placed randomly at the left x < 0 end of the system, see figure 1, the upper panel; periodic boundary conditions are used for the top and bottom edges. We investigate how cells move and proliferate to fill the void in the positive x direction.
In addition to random motion, a cell can perform directed motion along the gradient of a chemical; this phenomenon is called chemotaxis. We separately take into account positive chemotaxis, chemotaxis due to a chemoattractant (a cell moves toward the higher concentration of chemoattractant, e.g. nutrients) and negative chemotaxis due to a chemorepellent (a cell moves away from the higher concentration of chemorepellents, e.g. waste products). In the case of positive chemotaxis, the chemoattractant follows a standard diffusion equation with an extra consumption term where u is nutrient concentration, c is the consumption coefficient, D is the diffusion constant of chemoattractant and n is the cell density. Initially, we set a scaled nutrient concentration of u = 1 everywhere. We write down the simplest form of the consumption term: it is proportional to both cell density (more cells consume higher amounts of nutrients) and nutrient concentration.
No-flux boundary conditions are incorporated. For negative chemotaxis, the chemical is a waste product, which has a production rate c 1 and a breakdown rate c 2 : The chemorepellent production term depends only on the local cell density, and the usual linear degradation term (waste breakdown) is assumed. Notice that equation (2) can also be considered for positive chemotaxis models where the chemoattractant is generated by cells [23,24]. In the model, a cell is chosen randomly and tries to perform one of the three processes: proliferation, random hopping and directed migration (chemotaxis). As a result, a cell either proliferates or migrates into a neighboring lattice site. The probability of proliferation is given as a constant α. If the cell tries to hop (either performing random walk or directed migration), there is a possibility that adhesion to a neighboring cell (or cells) will prevent movement. If the adhesion does not prevent diffusion, then the cell will next decide which direction to jump. In the case of no chemotaxis, each direction simply has a 0.25 chance of being selected (if the corresponding neighboring site is not occupied). With chemotaxis added in, there is still a 0.5 chance that the cell will choose to jump along y versus jumping along x. The probability of jumping in +x direction versus −x direction (or +y versus −y) will be modified by the direction and magnitude of the chemical gradient at the cell's location. A positive chemotaxis coefficient makes the cells more likely to move toward high nutrient concentrations, while a negative value makes cells apt to move away. We assume that the cell can only sense the chemical concentration in the sites immediately surrounding it, and we calculate the effect of chemotaxis on the probability of hopping using the following: Here u is the chemical concentration (normalized to be between 0 and 1) and χ is the chemotaxis coefficient, given by χ = 2 c f /(1 + 3u) 2 , where c f is a (constant) chemotaxis factor. A positive value of µ x (µ y ) means that there is an increased probability of the cell jumping right (+x) (up, +y), while a negative value increases the probability of jumping left (−x) (down, −y). Therefore, for positive chemotaxis, c f > 0, so that cells are more inclined to jump toward the higher nutrient concentration. Likewise, for negative chemotaxis, c f must be negative. When χ (or δu) is zero, chemotaxis has no effect on cell movement. We consider the chemotaxis to be along the x-coordinate and denote µ = µ x . In our model, the overall probability of hopping is given by where q is the adhesion parameter and n is the number of neighboring cells (which can vary from 0 to 4 for a square lattice) [5,21]. Indeed, for high q, it is quite difficult for a cell to detach from its neighbors and migrate away. The migration includes random motion (diffusion) and chemotaxis. The probability of random hopping is then P diff = (1 − α) (1 − q) n /(1 + |µ|) (and then a cell jumps with probability of 1/4 to one of the neighboring sites), while the probability of directed hopping is P chem = (1 − α) (1 − q) n |µ|/(1 + |µ|) (and then a cell jumps with probability 1 along the chemical gradient). Thus, for |µ| 1 cells mostly perform random motion, while for |µ| 1, cells mostly perform directed motion. In simulations, |µ| 1 (therefore the term (1 + |µ|) can be omitted), but the effect of chemotaxis on front propagation is still significant, see below.
The expression χ = 2 c f /(1 + 3u) 2 [22] was found to qualitatively match the behavior of Escherichia coli bacteria; we assume that the same qualitative behavior is valid for glioma cells. Firstly, a larger chemical gradient δu produces a stronger effect. Secondly, as the chemical concentration increases, more and more chemoreceptors become occupied and eventually the cell will be unable to determine the direction of the chemical gradient. This explains an inverse relationship between the chemotaxis coefficient and nutrient density. This behavior is consistent with recent experimental results [25]. Researchers calculated the signal-to-noise ratio (SNR) determined by intracellular and extracellular stochastic fluctuations and showed that this parameter determines cell velocity due to chemotaxis [25]. It is shown in [25, figure 1] that higher values of SNR correspond to higher cell velocities (stronger chemotaxis), while [25, figure 3] shows how the SNR parameter depends on the gradient of chemical concentration and on the absolute value of chemical concentration. When the ratio of δu to u is constant, the cell velocity reaches maximum at some intermediate values of chemical concentration u, which is consistent with the form of chemotaxis coefficient introduced above.

Comparison with experiments
We performed simulations of the discrete lattice model and compared the results with the experimental observation on cell migration during 72 h. In figure 1, we see the snapshots of the system at the beginning and at the end of a typical simulation. Cells are initially placed at the left end. Over the course of the simulation, cells migrate and divide such that the area initially occupied has a much higher cell density and cells have moved right into the unoccupied area.
There is a significant discrepancy in the literature about the values of diffusion and degradation/production rates. Still, we have chosen a reasonable set of parameters to simulate chemotaxis. The diffusion coefficient of glucose in the brain is given by D = 6.7 × 10 −7 cm 2 s −1 = 67 µm 2 s −1 [26]. A cell jumps randomly its own diameter (20 µm) in time t diff (3 min). Therefore, cell diffusion coefficient in physical units is 400 µm 2 per 180 s. Therefore, the ratio of the two diffusion coefficients (or the scaled nutrients diffusion coefficient) is 30; this is exactly the diffusion constant for chemoattractants and chemorepellents used in simulations. The typical times for production and degradation were assumed to be of the order of 12 h (or a day). The corresponding scaled rates become t diff /12 h = 0.0042. This is the order of magnitude of the rates used in simulations. Figure 2 shows the chemical concentration as a function of the x coordinate along the system at the end of simulations. In the case of positive chemotaxis, the nutrient concentration was initially set to 1 everywhere. On the left end of the system, the concentration significantly decreased as the nutrients were consumed by the cells (which were initially placed at the left end of the system). Note that the furthest cell did not reach even 2000 µm, so the drop in nutrient concentration for most of the graph is due to diffusion of nutrients from the area without cells to the area with them. For negative chemotaxis, we have just the opposite. Initially, the chemical concentration was 0 everywhere, and is still near zero at the right end of the system. At the left end, cells produce waste, which then diffuses out to distances cells have not yet reached.
In figure 3, we compare the experimental observations [20] and the simulations of the lattice model both with and without chemotaxis. In the experiments, there is initially a mass of cells which is allowed to migrate/proliferate to the right for 72 h. Symbols show the results of three identical experiments. The snapshots of the system at different time points were analyzed by counting cells within the bins 110 µm wide. The simulations were set up to match the experiments as closely as possible, including the size of the system, the initial concentration of cells, the size of the bins in which the cells are counted, and the time given for activity. The simulations are averaged over ten runs to produce smoother curves and error bars.  (2) (for negative chemotaxis). The concentration was initially 1 everywhere for positive chemotaxis, and 0 for negative. This figure shows the average over the lateral coordinate y, and over ten runs. For positive chemotaxis, the consumption rate c was 0.0030, the chemotaxis factor was c f = 50, and the diffusion constant D was 30. For negative chemotaxis; the waste production rate c 1 was 0.0030, the waste breakdown rate c 2 was 0.0015, the chemotaxis factor was c f = 10 and the diffusion constant D was again 30. In the previous work [20], we looked at the same model without chemotaxis. We found out that the model without cell chemotaxis systematically underestimated the migration of the center of mass of migrating cells. Indeed, when chemotaxis is not included in the simulation, the cell density drops too fast compared to experiment, and the migration of the main mass of cells is underestimated, figure 3 (red solid line). Once we take into account chemotaxis, however, the simulations of the hybrid model produce cell density profiles consistent with the experimental results; this is true regardless of whether the migrating cells are being influenced by a chemoattractant (the results are shown by the blue dotted line) or chemorepellent (the results are shown by the green dashed line). The intuitive explanation is simple. For intermediate consumption (or production) rates, the location of the steepest gradient of chemoattractants (chemorepellents) corresponds to intermediate density of cell, so there the effect of chemotaxis is the largest. In contrast, cells in the leading edge of the front 'see' shallow chemical gradients and therefore their migration is almost unaffected by chemotaxis.
Positive and negative chemotaxis provide similar results since both phenomena lead to increased migration along the gradient of some chemicals. This means that one cannot distinguish between the two mechanisms focusing just on cell density profiles without knowing the nature of chemoattractants or chemorepellents. It has been found that large Slit proteins secreted by cells may function as chemorepellents binding Robo receptors [27]. To investigate the role of Slit2 in cell migration, we propose an experiment where a cluster of cells is placed  on a substrate with a nearby chemical source of Slit2. One can divide the substrate into four quadrants, with the chemical source in just one of those quadrants. By comparing the migration of cells in that quadrant to the others, we can distinguish between directed cell migration and random cell motion. This type of experiment has been done before to show that stromal cell-derived factor 1α is a chemoattractant for neural progenitor cells (see [28, figure 2]). We also suggest to vary the concentration of Slit2 to verify the saturation effect: when a high concentration of a chemical is introduced, it may bind with all available Robo receptors on the cells, making it impossible to detect the direction of the chemical gradient.
Once it is clear that Slit2 is a real candidate for a chemorepellent factor, one can try to distinguish the positive and negative chemotaxis mechanisms. We propose to place a tumor spheroid on a substrate, where the gradient of cell culture media is established. Negative chemotaxis mechanism will lead to spherically symmetric invasion; the deviations from spherical symmetry will show the relative role of positive chemotaxis.   (1)) and the continuous model (equations (1) and (6)). Two different times (t 1 = 6000 cell diffusion times and t 2 = 12 000) are shown, the coordinate x is measured in cell diameters. The cell density, n, is shown as a solid red line for the simulations of the continuous model and blue circles correspond to the hybrid simulations. The corresponding nutrient profiles are shown as a dashed and dotted lines, respectively. For these simulations, the adhesion parameter q = 0, nutrient consumption rate c = 0.06 and the chemotaxis factor c f = 100.

Derivation of the macroscopic continuum model
In order to analyze the phenomenon of front propagation theoretically, we derive the continuum equation for cell density from the underlying microscopic lattice model. Firstly, we write down the phenomenological master equation for P(i, j), the probability that the site i, j is occupied at time t. Schematically, the equation is written aṡ P(i, j) = hopping due to diffusion + hopping due to chemotaxis + proliferation.
The terms 'hopping due to diffusion' and 'hopping due to chemotaxis' in equation (4) consist of several terms describing jumps into the site (i, j) and out of that site. The diffusion and proliferation terms in the case of nonzero adhesion are derived in [20,29]. Assuming chemotaxis in the positive x-direction, we get the following chemotaxis terms:  (1)) and the continuous model (equations (1) and (6) The first term describes a jump from site i − 1, j to site i, j. For this event to occur, the site i − 1, j should be occupied and the site i, j should be empty. Note that the joint probability of these two events is written as a product P(i − 1, j)(1 − P(i, j)), since we neglect correlations between P(i − 1, j) and P(i, j); this assumption is valid only for low adhesion (small q) as well as for small proliferation. Note that for large adhesion the probabilities that the two neighboring sites are occupied are not independent; in particular, supercritical adhesion leads to clustering [5]. Correlations cannot be neglected also when the proliferation rate is much larger than the diffusion rate, this case leads to a completely different dynamics and the velocity of front propagation differs significantly from the one predicted by the Fisher-Kolmogorov equation [30]. We have also assumed that |µ| 1, although the derivation can be easily modified for the general case.
Taylor-expanding these probabilities up to the first order terms, we arrive at the following chemotaxis term, which also accounts for cell-cell adhesion: where n is a cell density and the chemotaxis coefficient is given by µ = [2 c f /(1 + 3u) 2 ] (∂u/∂ x). The full equation for cell density becomes (compare with [20, equation 4]): where D(n) = (1 − q) 4n [1 + 4n (1 − n) ln(1 − q)]/4 [29], the time is measured in the units of cell diffusion time t diff , and the distance is measured in the units of cell diameter.

Front propagation: adhesion and chemotaxis
Now we employ both the discrete lattice model and the derived continuum equation (6) to investigate the phenomenon of front propagation. The lattice model was simulated as described in section 2; to obtain the cell density profile as a function of x-coordinate, we performed averaging over the lateral direction and over many runs. In order to identify the concentration of chemicals around a cell in the lattice model, the equation for chemoattractants (or chemorepellents) were solved numerically using the third-order Runge-Kutta method. Within the theoretical continuum framework, equations (6) and (1) or (2) were solved using the Matlab package. Figure 4 shows front propagation in the case of a positive chemotaxis: cells consume nutrients and (in addition to random motion) they migrate toward the higher concentration of nutrients. The propagating fronts (both of cells, n(x) and of nutrients, u(x)) are computed both by employing hybrid discrete-continuum approach and by numerically solving equations (1) and (6). The agreement is very good. Notice the step-like density profile that can be observed for high enough chemotaxis coefficient (or for high enough consumption rate). This feature may represent discontinuity in density gradient; a similar observation was done already in a simpler setting without adhesion and without diffusion of nutrients for high enough chemotaxis coefficient [31]. The difference in cell density profiles is clearly seen in figure 5, where the front profile n(x − vt) is plotted in cases (q = 0 and µ = 0), (q = 0 and µ = 0) and a general case (q = 0 and µ = 0). Interestingly, nonzero adhesion seems to substantially smooth the density step. Notice an excellent agreement between the simulations of the discrete lattice model (symbols) and the continuum approach (lines). Figure 6 shows the velocity of front propagation, again computed both by simulating stochastic lattice model for cell dynamics (symbols) and by solving equations (1) and (6) numerically (lines). One can see an excellent agreement for zero adhesion and some discrepancy for nonzero adhesion. We believe that this discrepancy occurs due to neglecting correlations in the derivation of equation (6) from the microscopic lattice model. The speed of cell invasion increases with the consumption rate c since the nutrients profile becomes sharper for larger c, increasing the effect of chemotaxis.

Conclusions
In this paper we studied the effect of chemotaxis on migration of proliferative and adhesive cells. We have derived a macroscopic continuum equation for cell density that accounts both for chemotaxis and cell-cell adhesion. Numerical solution of this equation (coupled to the equation for chemical concentration) shows a good agreement with the simulations of the underlying stochastic lattice model. Both the front profiles and the front velocity are correctly captured by the derived equation. It is expected that some discrepancy should occur for intermediate values of adhesion parameter, since it introduced correlations that were not taken into account in the continuum approach. One way to incorporate these correlations is by using a moment closure approach [32]. Figure 5 shows a step-like cell density profile when the chemotaxis is strong enough. The exact nature of this 'step' will be considered in the future. While clearly still there, the addition of cell adhesion to the models caused a substantial decline in the prominence of this step. This is intuitively clear since adhesion decreases the effect of chemotaxis and the effective chemotaxis coefficient in equation (6).
We have also shown that both positive chemotaxis (where cells are driven toward a chemoattractant, e.g. a nutrient) and negative chemotaxis (where cells are driven away from a chemorepellent, e.g. a waste product) can explain recent experimental observations of migration of glioma cells on a substrate [20]. The addition of chemotaxis makes the main mass of cells migrate faster, which widens the front profile and causes the agreement with experiment that was lacking in earlier models. Since there was no external gradient of growth factors and there was a constant supply of nutrients in experiments [20], it seems probable that cells secreted waste chemicals and then tried to move away from that region similar to what was observed in [15], but the exact chemorepellent factor remains to be found. One possible chemorepellent factor that was found to be important in glioma cell migration is a Slit protein that binds to Robo receptors expressed by glioma cells [27].