How cells stay together; a mechanism for maintenance of a robust cluster explored by local and nonlocal continuum models

Migration of clusters of tumor cells found by [16] introduced the importance of collective cell behaviour and migration. However, cell clusters form (and migrate collectively) in many normal developmental and morphogenetic processes [37]. This is true of the post-aggregation “slug” stage of social amoebae such as Dichtyostelium discoideum (DD) [7] and the primary lateral line primordial (PLLP), a cluster of about 100 cells that deposits sensory organs along the length of a zebrafish embryo [9]. Long range migration by a cluster of cells is essential in the formation of the sympathetic nervous system, as studied in [20] for chick embryos. Cell clusters can depend on actual contact, adhesion, and short or long-lived junctions formed between cells [15]. But in the case of the primary sympathetic ganglion (PSG) cluster described in [20], the cells form a “swarm” and avoid adhering to one another.

A question arising in these experiments is what determines the observed cell cluster sizes, and in particular how cluster size is related to cell-cell interactions. From spatial RNAseq experiments in the PSG, it is known that the cells in the migrating cluster interact via chemical signalling [20], but whether the chemicals are attractants, repellents, or both, is unclear. More generally, an interesting biological problem is to determine what combination of such chemicals (with various rates of diffusion, secretion, and decay) would be consistent with the formation and maintenance of robust cell clusters. This general question motivates our paper.

Models of collective cell behaviour have utilized the cellular Potts Model [29] and a spheroidal (3D) cell-based computation of the DD slug [35], and of the PLLP [23], and other custom-built 2D cell migration models for neural crest cells [27]. See review in [4]. Comparison of continuum and agent-based models in cell populations are found in [3, 11].

Non-local equations have frequently been used to describe cell densities. These include models that have non-local interactions a-priori such as the cell-cell adhesion model of [1], the non-local chemotaxis model by [32], or models in which the non-local term appears due to a “model reduction” technique (e.g., quasi steady-state assumptions, as in [23, 28]). A popular reduction of non-local models to local models has been considered in the literature by [30], by [17], and recently by [14]. In this reduction, the limits of integration are formally taken to zero, so that fully local equations are found. These local equations can give powerful insights and connections to well studied local equations e.g. Cahn-Hilliard or Allen-Cahn [8].

While much of our paper rests on the results of [14], we here pose several new questions of both mathematical and biological flavour:

1. (a)

   What biological processes give rise to cell-cell interactions? (Rephrased, where do the kernels describing cell “potentials” come from?)

2. (b)

   Under what conditions (on the kernels or the parameters) would a single cell species create a stable cluster that is maintained?

3. (c)

   How does the cluster size (radius) depend on the cell interaction parameters?

4. (d)

   For a given underlying cell signaling model, can we create a “map” in parameter space that shows the possible regimes of behaviour of a population of such cells? What behaviours can occur other than clustering?

5. (e)

   How does a local (PDE) approximation compare with the full nonlocal (integro-PDE) model?

6. (f)

   How do the predictions of the continuum models compare with agent-based models with similar mechanisms?

More broadly, we espouse the idea of using several kinds of model simplifications to gain initial understanding of processes, before attempting more complex scenarios and adding details. Indeed, we make several simplifications to obtain insight on the migration of the PSG cluster: (1) We first drop the analysis of migration and consider the problem of cluster existence, shape, and stability. (2) We carry out analysis in a 1D geometry, where results are simplest and most intuitive. (3) We assume one species of identical cells.

In what follows, we will examine both nonlocal models and their local approximation, applying previous results of [14, 28] to the questions posed above.

To keep this paper self-contained, we first provide a pedagogical summary of part of [14], listing some basic results from that paper.

Let $\rho(x,t)$ denote density of cells, at position $x$ and time $t$. Ignoring random cell motion, the full nonlocal model in generality is

where $W$ represents a potential induced by a cell on its neighboring cells at some position $\mathbf{x}$ at time $t$. The velocity $\mathbf{v}$ is proportional to a force $\mathbf{F}$ in the low Reynold’s number regime of cells; here the drag coefficient has been scaled to 1. The force is assumed to be the gradient of a potential induced by the cells. $K*\rho$ is a convolution that sums up all distance-dependent cellular interactions.

Here we consider the spatial 1D case, where the model reduces to

where

is a convolution over the whole real line, $-\infty<x<\infty$. We will take the simplest case, where $K$ is a symmetric kernel, i.e. cells interact in a symmetric way, making $K$ an even function of space, $K(x)=K(-x)$, so cells induce a potential that is symmetric about their position. This general problem was one of the topics studied in [14]. We wrote the convolution in two forms, as the latter is useful to our next step.

As in [14], we consider a local approximation for the nonlocal 1D model (2). To do so, we transformed the integration variable ($z=x-s$) as shown above. Expand the term $\rho(x-z,t)$ in a Taylor series about $x$ to obtain a sum of partial derivatives $\rho,\rho_{x},\rho_{xx}$ etc. The coefficients of such terms are proportional to moments of the kernel, that is

The even symmetry of the potential implies that $a_{1}=0$. We keep terms up to order $z^{2}$. Then the local approximation is

We apply the above general theory of [14] to a biologically relevant setting where cell signaling is based on known or plausible mechanism(s). We here introduce the Morse potentials, and then show that they capture known cell-cell signaling mechanism(s).

The Morse potential we consider combines repulsion and attraction, written as

where $A,R\geq 0$ are (attraction, repulsion) force magnitudes and $a,r\geq 0$ are typical spatial scales of attraction and repulsion. This is an even function. It is to be interpreted as a symmetric “potential” induced by a single cell located at the origin on its neighbors at another location $x$. We have used the above notation to be consistent with [28] where forces (rather than potentials) were employed in a similar setting. It will be convenient to define the following dimensionless ratios:

Our results will be summarized as regimes in the $C\ell$ plane, where inequalities determine certain regions whose boundaries are expressed in terms of $C$ and $\ell$.

The Morse potential was introduced in 1926 by Morse to replace the harmonic potential in quantum mechanics systems for diatomic molecules. The Morse potential proved to be more accurate given experimental observations. Interestingly no derivation of this potential from first principles was provided, to our knowledge, until [28], and much later in [10]. Below, we briefly demonstrate the chemical attraction-repulsion system of [28] that directly leads to a Morse type potential.

Following [28], we consider a direct generalization of the model of [22] that includes chemotaxis up the gradient of a diffusible attractant $c_{a}(x,t)$, and down the gradient of a diffusible repellent, $c_{r}(x,t)$, both secreted by the cells. The PDEs describing the cell density $\rho(x,t)$ and the chemical concentrations are

Here $D_{\rho},D_{a},D_{r}$ are, respectively the cell random motility and the chemical rates of diffusion. The parameters $\chi_{a},\chi_{r}$ are chemotactic coefficients for cells moving towards (away from) the attractant (repellent). The parameters $s_{a},s_{r}$ are rates of secretion of the two chemicals by cells and $k_{a},k_{r}$ are rates of chemical decay. The parameters $\epsilon_{a,r}$ are assumed to be small, so that the chemical dynamics operates on a faster timescale than the cell motion.

Signaling molecules diffuse at a far faster rate than the motion of cells. Hence, a quasi steady state (QSS) approximation for the two chemical PDEs can be justified by time-scale separation, implying that

Eqn (24) can be solved by the Green’s function method [28], leading to

We recognize the above as a convolution of the cell density with a kernel $K_{j}(x)$ of the form

Here the expression $p_{j}$ is the spatial range of the kernel $K_{j}(x)$ (corresponding to $j=a,r$ for the attraction and repulsion). In our case, there is a superposition of two such terms, one each for the attractant and the repellent. The quantity $P_{j}$ will prove to be proportional to the amplitudes of the kernel ($A,R$), but we include the chemotactic sensitivities of the cells in those amplitudes, see below.

We use these QSS solutions to eliminate $c_{a},c_{r}$ from the PDE system, but first we rewrite the cell PDE in the compressed form

We chose signs above so that $W$ can be interpreted as a chemical potential; for example, if there is no repellent ($c_{r}=0$) then there would be a “potential well” where the attractant is most concentrated.

The form of the expressions in (26) results in spatial ranges of the repellent and attractant, and attraction repulsion force magnitudes

Then the chemical potential can be written as

The above implies that $W$ acts as a net “chemical potential landscape” induced by the cell density due to chemical signaling. This potential sets up the attraction-repulsion force-field that results in cell motion. (As the cells move, that potential landscape will shift.)

Substituting the values of $r,a,R,A$ from (28) into the conditions (22) for a stable cluster, we find that

These can be combined into a single statement, namely that for a robust cluster to form,

The inequality (30) is the only requirement for the existence of a (stable) cluster.

If we normalize all parameters by the rates for the attractant, then we find that the relative diffusion of the repellent should be smaller than the product of cell’s repulsion taxis rate and its repellent secretion rate. The repellent relative decay rate should be larger than both of the above. In other words, for a cell cluster to form, the repellent should diffuse slowly and decay rapidly so that its spatial range is small; the cell should have a repulsion taxis and repellent secretion rate that are neither too large nor too small relative to the attractant.

We can see from the above that the magnitudes of repulsion/attraction forces $R,A$ can be “controlled” by rates of secretion $s_{j}$ of the repellent and attractant as well as the cell’s chemotactic sensitivity to gradients of each chemical, i.e. the parameters $\chi_{j}$. The spatial ranges of the interactions, $r,a$ are set by the diffusivity and the decay rates of the chemicals. The lower the rate of diffusion, the smaller the spatial range, and the larger the magnitude of the respective attraction-repulsion force. The decay rate can be adjusted if cells also secrete enzymes that degrade the ligands. For example, the social amoebae Dictyostelium discoideum secrete an enzyme, phosphodiesterase, that degrades the attractant cAMP.

%ֿsectionOther models

While the above theory proves particularly “elegant” for Morse potentials, it has a straightforward generalization to other cases, where the repulsion is not due to chemical signaling. We discuss one example, in which the repulsion is due to a power-law force term, analogous to the concave down portion of a double-well potential. Suppose the repulsive potential is given by

where $r$ is the spatial range and $R$ the force magnitude, as before. This is a localized compact repulsion, akin to a non-rigid volume exclusion.

We can find the analogous PDE approximation as before, by computing the moments of the repulsive potential (31), to be

where some constants, ($\gamma_{0}=16/15,\gamma_{2}=8/105$) now appear from the polynomial integration step, as shown in the Appendix B. We observe that the results are similar to those of Morse kernels, save for some multiplicative constants.

Now combining this repulsion with the previously described chemoattractant attractive potential, we obtain similar features, slightly modified by the presence of two constants. For example, we find that the spatial frequencies, previously given by (17) are now given by

(The factors of 2, present in $a_{0},a_{2}$ for the double-Morse potential case in (13) no longer cancel out, and hence appear only in the attraction.) The dimensionless cluster size is $\mathcal{F}=\pi/\tilde{\omega}$. The analysis carries over, but the boundaries of the behavioural regions in the $C\ell$ plane are now defined by curves $C=2\ell^{2}/\gamma_{0}$, $C=2\ell^{4}/\gamma_{2}$, etc. We can see the effect of this distinct repulsive potential by considering the $C\ell$ plane on a log-log plot where the curves are of the form

where $G_{0}=\log(\gamma_{0}/2),G_{2}=\log(\gamma_{2}/2)$ are merely constants that shift the boundary curves up or down. (See right panel of Fig. 2.) The new repulsive potential also results in changing of the size of the cluster (where a cluster exists), shown by the somewhat distinct hues on that panel compared to the double-Morse potential case.

More generally, one can consider a repulsive potential of the form

for some generic function $\phi$. Provided the function $\phi$ is even, integrable on $-\infty<x<\infty$, the same general conclusions hold, (given that $\phi(x)$ is normalized) but with somewhat specific constants $\gamma_{0},\gamma_{2}$ that result in shifting boundary curves up or down in the $C\ell$ plane. In the above example, the repulsive potential has compact support $-r<x<r$, but this need not be the case in general, as we have seen with Morse repulsion.

So far, the local and nonlocal Morse model analysis informs us about conditions for instability of a uniform density solution (HSS) and about conditions for existence of a finite cluster. However, we still lack two important pieces of information: (1) Conditions that guarantee a stable cluster and (2) an assessment of the accuracy of the local approximation of the full nonlocal problem.

We sought to address these issues by numerically solving both the local and the nonlocal models in 1D, with the Morse potential specified in Eqn. (12). For the non-local model, we used a domain far larger than the expected cluster size and periodic boundary conditions. We simulate the problem

with $K$ the Morse potential whose kernel is given by (12). For details on the numerical methods see Appendix A.

In this paper, we have addressed several questions. First, we examined one example of cellular signaling via diffusible chemotactic chemicals that leads to effective nonlocal cell-cell interactions. This example was first proposed by [28], but has not yet been recognized in recent nonlocal modeling efforts. We showed that when cells secrete a diffusible chemical that decays, and that attracts or repels the cells via the process of chemotaxis, it essentially sets up a “chemical potential” of the Morse type.

Morse potentials have been used in more general settings of cell-cell mechanical interactions in computational biology papers. Morse potentials appear in the context of cell adhesion models [1], as briefly surveyed in Appendix  C. [25] used a Morse potential to model adhesive/repulsive interactions of red blood cells (RBCs) in an advanced computational model of RBC aggregation using an immersed finite element method. [39] used a modified Morse potential to represent adhesion of platelets to one another in an ABM model for blood clotting. The Morse potential was adapted to experimental data on shear conditions. [21] used cell tracking in mouse embryos and C. elegans to infer effective mechanical forces and pairwise potential of cell-cell interactions. They showed that Morse potentials were generally a good fit, since tuning the Morse parameters allowed for a range of interactions. They also showed that species differences that implied distinct potential profiles also accounted for tight vs loose cell aggregate morphologies. As discussed in [31], many forms of potentials are reasonable, provided they account for four parameters for magnitudes and spatial ranges of attraction and repulsion, akin to our $A,R,a,r$. Hence, even without the direct link to mechanism, Morse potentials have proven to be convenient approximations for a variety of non-local interactions.

A second question that we addressed is what kinds of attributes chemotactic signals should have to create a robust cell cluster. We showed that the behaviour of cells with attractant and repelent interactions can be classified in a $C\ell$ parameter plane where these dimensionless parameters are ratios of magnitudes and spatial ranges of the attractant and repellent. Similar curves were found in [28] in connection with pattern formation criteria. We extended that theory to predict the size of the stable cluster in terms of the chemical properties. We also found the variety of possible behaviours for regions other than that of the cluster formation. In most such regions, the cells cannot aggregate due to dominance of repulsion or the cells collapse due to dominance of attraction.

Some of the theory and predictions here and in [14] rely on the ability to approximate a nonlocal problem by a local PDE, where steady state solutions can be written explicitly and conveniently (at least in 1D). We quantified the degree of accuracy of such local approximations by numerically simulating and comparing the predictions of the nonlocal problem and its local PDE approximation.

Another question of interest is how the continuum results compare with agent based models. Analytic results for the case of a well-spaced swarm (where all agents are at equal distances from neighbors) in [26] arrived at similar regime boundaries in the $C\ell$ plane, but did not explore the regime of a loose cluster of finite radius. In [12], the authors investigated an agent-based model for velocity-limited particles that interact with Morse attraction-repulsion kernels. They classified behaviours into “H-stable” (well-spaced group) and “catastrophic” (collapse into a tight cluster), on a parameter plane very similar to the $C\ell$ plane shown here. The authors of [24] consider asymptotic dynamics of a 1D swarming model with attraction-repulsion kernels, including the Morse kernel. They show how properties of the kernel’s moments account for whether the group will contract or spread, or form a steady-state cluster. In [36], the authors include self-propulsion, friction, together with attractive/repulsive Morse potential. They assume that the total mass stays constant, regardless of the size of the flock. They numerically simulate their model in 3D, and find regimes for clumps spheres, dispersion, mills, rigid-body rotation, flocks. The above examples have features that go beyond our continuum models, and yet some of the properties that we find persist in these distinct situations.

In future plans for the analysis of collective behaviour, and particularly for the migration of the presypathetic ganglion (PSG) of [20], it would be important to consider the interactions of a second cell type (as modeled in [14]) with similar signaling mechanism, as well as the migration of a cluster towards its target. In such extensions, some analysis is available, though results are not as easily stated, as many new parameters enter the models.