Jeans type instability for a chemotactic model of cellular aggregation
We consider an inertial model of chemotactic aggregation generalizing the Keller-Segel model and we study the linear dynamical stability of an infinite and homogeneous distribution of cells (bacteria, amoebae, endothelial cells,…) when inertial effects are accounted for. These inertial terms model cells directional persistance. We determine the condition of instability and the growth rate of the perturbation as a function of the cell density and the wavelength of the perturbation. We discuss the differences between overdamped (Keller-Segel) and inertial models. Finally, we show the analogy between the instability criterion for biological populations and the Jeans instability criterion in astrophysics.
pacs:05.45.-aNonlinear dynamics and nonlinear dynamical systems
The self-organization of biological cells (bacteria, amoebae, endothelial cells,…) or even insects (like ants) due to the long-range attraction of a chemical (pheromone, smell, food,…) produced by the organisms themselves is a long-standing problem in physical sciences murray (). This process is called chemotaxis. The chemotactic aggregation of biological populations is usually studied in terms of the Keller-Segel model ks ():
which consists in a drift-diffusion equation (1) governing the evolution of the density of cells coupled to a diffusion equation (2) involving terms of source and degradation for the secreted chemical . The chemical is produced by the organisms (cells) at a rate and is degraded at a rate . It also diffuses according to Fick’s law with a diffusion coefficient . The concentration of cells changes as a result of an oriented chemotactic motion in a direction of a positive gradient of the chemical and a random motion analogous to diffusion. In Eq. (1), is the diffusion coefficient of the cells and is a measure of the strength of the influence of the chemical gradient on the flow of cells. These coefficients depend a priori on the concentration of cells and on the concentration of the chemical. The Keller-Segel model is able to reproduce the chemotactic aggregation (collapse) of biological populations when the attractive drift term overcomes the diffusive term .
However, recent experiments of in vitro formation of blood vessels show that cells randomly spread on a gel matrix autonomously organize to form a connected vascular network that is interpreted as the beginning of a vasculature gamba (). This phenomenon is responsible of angiogenesis, a major actor for the growth of tumors. These networks cannot be explained by the parabolic model (1)-(2) that leads to pointwise blow-up. However, they can be recovered by hyperbolic models that lead to the formation of networks patterns that are in good agreement with experimental results. These models take into account inertial effects and they have the form of hydrodynamic equations gamba ():
The inertial term models cells directional persistance and the general density dependent pressure term can take into account the fact that the cells do not interpenetrate. In these models, the particles concentrate on lines or filaments gamba (); filbet (). These structures share some analogies with the formation of ants’ networks (due to the attraction of a pheromonal substance) and with the large-scale structures in the universe that are described by similar hydrodynamic (hyperbolic) equations, the Euler-Poisson system peebles (). The similarities between the networks observed in astrophysics (see Figs 10-11 of vergassola ()) and biology (see Figs 1-2 of gamba ()) are striking.
including a friction force . This type of damped hydrodynamic equations was introduced in csr (); gt () at a general level. This inertial model takes into account the fact that the particles do not respond immediately to the chemotactic drift but that they have the tendency to continue in a given direction on their own. However, after a relaxation time of order , their velocity will be aligned with the chemotactic gradient. This is modeled by an effective friction force in Eq. (7) where the friction coefficient is interpreted as the inverse of the relaxation time. This term can also represent a physical friction of the organisms against a fixed matrigel. In the strong friction limit , or for large times , one can formally neglect the inertial term in Eq. (7) and obtain gt (); crrs (); virial ():
Substituting this relation in Eq. (6), we recover the parabolic Keller-Segel model (1)-(2). Therefore, the Keller-Segel model can be viewed as an overdamped limit of a hydrodynamic model involving a friction force. Alternatively, neglecting the friction force , we recover the hydrodynamic model introduced in gamba ().
In this paper, we study the linear dynamical stability of an infinite and homogeneous distribution of cells with respect to the inertial model (6)-(8). We determine the condition of instability and the growth rate of the perturbation, and discuss the differences with the results obtained with the Keller-Segel model (1)-(2). We also discuss some analogies with the dynamical stability of self-gravitating systems. Indeed, there are many analogies between the chemotactic aggregation of biological populations and the dynamics of self-gravitating Brownian particles crrs (). In particular, the Keller-Segel model (1)-(2) is similar to the Smoluchowski-Poisson system crs () and the hydrodynamic equations (6)-(8) are similar to the damped Euler equations of Brownian particles gt (); virial (). The main difference between biological systems and self-gravitating Brownian particles is that the Poisson equation in the gravitational problem is replaced by a more general field equation (8) taking into account the specificities of the biological problem. Owing to this analogy, we shall discuss the relation between the instability criterion of biological populations and the Jeans instability criterion jeans () in astrophysics.
2 Instability criterion for biological populations
2.1 The dispersion relation
Linearizing the equations around this stationary solution, we get
Looking for solutions of the form and , we get
where . These equations have non-trivial solutions only if the determinant of the system is equal to zero yielding the dispersion relation
The condition of marginal stability () corresponds to and the condition of instability is
We note that the instability criterion does not depend on the value of and . Let us consider in detail some particular cases.
2.2 The case
The discriminant is
and the two roots are
If the perturbation decays exponentially rapidly and if the perturbation grows exponentially rapidly. In that case, the system is unstable and is the growth rate of the perturbation. If , then so that (unstable). Alternatively, if , then . Either and (stable) or and (stable). Therefore, the system is unstable if
and stable otherwise. To determine the range of unstable wavelengths, we must consider different cases:
2.2.1 If :
In that case, a necessary condition of instability is that
If this condition is fulfilled, the unstable wavenumbers are determined by
The growth rate of the perturbation with wavenumber is . Therefore, the most unstable mode is the one which maximizes . It is given by
The largest growth rate is then determined by
In particular, for , we have
The range of unstable wavelengths is determined graphically in Fig. 1 and the evolution of the growth rate of the perturbation as a function of the wavenumber is plotted in Fig. 2 (we have also consider the case in this Figure by solving Eq. (2.1) which is a second degree equation in ).
2.2.2 If :
In that case, . The unstable wavenumbers are determined by
The most unstable mode is and the largest growth rate is given by
2.2.3 If :
In that case, the system is unstable for
The growth rate diverges when
corresponding to . Close to the critical wavenumber , we have
This expression is valid for finite and . It can be directly obtain from Eq. (19) by using for . Thus, when the temporal term is neglected in Eq. (8), i.e. , a critical behaviour occurs. This critical behaviour is regularized for . Indeed, taking , i.e. , in Eq. (2.1), we obtain
Taking the limit , we find that
which is finite for but diverges like when . For and , the dispersion relation can be simplified in
On the other hand, for , Eq. (2.1) reduces to
The positive root of this equation is
which is finite for but diverges like when .
The range of unstable wavelengths is determined graphically in Fig. 3 and the evolution of the growth rate of the perturbation as a function of the wavenumber is plotted in Fig. 4 (we have also consider the case in this Figure by solving Eq. (2.1)). We see that the case is very special. For , the range of wavenumbers seems to be stable according to Eq. (30) because the two roots of Eq. (19) are negative. However, for any finite value of , a third root appears. This root is positive and tends to infinity when (see Eqs. (34) and (38)). Therefore, this unstable branch is rejected to infinity when . This implies that the region is in fact extremely unstable for . In particular, for , the most unstable mode is and the largest growth rate is given by Eq. (38).
2.3 The case
In the overdamped limit, the hydrodynamical equations (6)-(8) return the Keller-Segel model (1)-(2). Let us consider the stability analysis in that case for comparison with the inertial case. The dispersion relation now reads
and the two roots are
Writing the solution in the form
it is easy to see that the system is stable if (i) () and if (ii) (). It is unstable otherwise. Since (ii) implies (i), the system is stable if and unstable otherwise. Thus, the system is unstable if
and stable otherwise. The range of unstable wavelengths is determined graphically in Fig. 5.
2.3.1 If :
In that case, the system is unstable for
and stable otherwise. This is the same criterion as for the inertial model. A necessary condition of instability is that
If this condition is fulfilled, the unstable wavenumbers are such that
These results are unchanged with respect to the inertial case.
We now determine the value of the optimal (most unstable) wavenumber and the corresponding growth rate (see Fig. 6). Eq. (39) is a second order equation in whose solutions are given by Eq. (40). We can maximize to obtain and . However, it appears simpler to proceed differently. Eq. (39) can also be viewed as a second order equation in of the form
There will be two roots and provided that and . This last condition can be written
The discriminant of Eq. (51) is given by
The condition to have two roots and is equivalent to with
(Note that the possibility must be rejected since it does not satisfy the requirement ). For , the two roots coincide so that is the maximum growth rate. It is reached for an optimal wavenumber , i.e.
2.3.2 If :
2.4 The case and
If we neglect the temporal term in the Keller-Segel model (), we obtain the dispersion relation
so that is explicitly given by
The instability criterion is given by Eq. (22).
2.4.1 If :
This is a particular case of Sec. 2.2.1 corresponding to . The expression of the largest growth rate is given by
The other results are unchanged.
2.4.2 If :
The system is unstable for the wavenumbers determined by Eq. (30). For , corresponding to , the growth rate diverges like
This divergence is regularized if . Taking , i.e. in Eq. (39), we get
For , we obtain
which is finite for but diverges like when . For and , Eq. (39) can be simplified in
On the other hand, for , Eq. (39) leads to
which is finite for but diverges like when .
Equations (32)-(34) and Eqs. (62)-(64) differ because both and tend to infinity, so that the expression depend on how the limits are taken. The general case can be treated as follows. For , taking in Eq. (19) we get
On the other hand, for , taking (i.e. ) in Eq. (2.1), we get
Finally, for and , we have
3 Analogy with the Jeans problem in astrophysics
3.1 The damped Euler equations
For , we can neglect the inertia of the particles so that . Substituting this relation in Eq. (71), we obtain a special case of the Keller-Segel model
where we have set . Equations (71)-(72) can be viewed as fluid equations appropriate to the chemotactic problem. Equation (71) is an equation of continuity and Eq. (72) is similar to the Euler equation where plays the role of a pressure and the chemotactic attraction plays the role of a force. Since , these equations describe a barotropic fluid. The main novelty of these equations with respect to usual hydrodynamical equations is the presence of a friction force which allows to make a connection between hyperbolic () and parabolic () models.
This hydrodynamic model including a friction force is similar to the damped barotropic Euler-Poisson system which describes a gas of self-gravitating Brownian particles gt (); virial () or the violent relaxation of collisionless stellar systems on the coarse-grained scale in astrophysics csr (). In that analogy, the concentration of the chemical plays the role of the gravitational potential . The main difference between the two models is that the Poisson equation for self-gravitating systems is replaced by a more general field equation (73) for bacterial populations. To emphasize the connection with astrophysical problems, let us consider a particular case of Eqs. (71)-(73) where , and and are constant. Then, introducing notations similar to those used in astrophysics (noting , , ), we can rewrite Eqs. (71)-(73) in the form
which can be interpreted as a generalized Smoluchowski equation. Thus, for , we obtain the generalized Smoluchowski-Poisson system describing self-gravitating Brownian particles in an overdamped limit gt (). Alternatively, for we recover the barotropic Euler-Poisson system that has been studied at length in astrophysics to determine the period of stellar pulsations cox () and the formation of large-scale structures in cosmology peebles (). Therefore, the chemotactic model (75)-(77) is similar to astrophysical models with additional terms. In the astrophysical context, the case would correspond to a shielding of the gravitational interaction on a typical length . This Yukawa shielding appears in theories where the graviton has a mass but in that case is very small which does not need to be the case in the biological problem.
We now consider the linear dynamical stability of an infinite and homogeneous solution of Eqs. (75)-(77). By mapping the equations (75)-(77) onto a generalized astrophysical model, we shall see that the instability criteria obtained in Sec. 2 are connected to (and extend) the Jeans instability criterion of astrophysics jeans (). To avoid the Jeans swindle bt () when , we have introduced a “neutralizing background” in Eq. (77). In fact, in the biological problem, this term appears naturally when we consider the limit of large diffusivity of the chemical (see jager ()); therefore, there is no “Jeans swindle” in the biological problem based on Eq. (2). A similar term appears in cosmology when we take into account the expansion of the universe and work in the comoving frame peebles ().
3.2 The dispersion relation
For a pure Newtonian interaction with , we have . Linearizing the equations around this stationary solution, we get
Looking for solutions of the form , we get
From these equations, we obtain the dispersion relation
In the case and , we recover the usual Jeans dispersion relation jeans ():
3.3 Instability criterion
If we set , the dispersion relation becomes
The two roots are
Accordingly, the system is unstable if
and stable otherwise. A necessary condition of instability is
If this condition is fulfilled the unstable wavelengths are such that
The wavelength which has the largest growth rate is given by
and the corresponding growth rate is given by
The instability criterion (91) is equivalent to the instability criterion (22) of Sec. 2.2 for the particular case of the chemotactic model considered in Sec. 3.1. The parallel with astrophysics is interesting to develop in order to give a more physical interpretation to Eq. (22). All the other formulae can be interpreted accordingly. In particular, we note that the coefficient plays the role of a velocity of sound .
3.4 Particular cases
Let us consider particular cases of the foregoing expressions:
For , one has , and
For (Newtonian potential), one has , (Jeans length), and
For , one has
For , one has
3.5 Isothermal gas
For an isothermal gas with an equation of state (for simplicity, we have noted instead of ), the velocity of sound is equal to the square root of the temperature: . It is relevant to re-express the previous relations as follows. Introducing the critical temperature
the growth rate of the perturbation can be written
The condition of instability reads
and a necessary condition of instability is . For the unstable wavenumbers (see Fig. 8) are such that with
The wavenumber with the largest growth rate is given by
and the largest growth rate by
The maximum wavenumber and the most unstable wavenumber are plotted as a function of the temperature in Fig. 9. In Fig. 10, we represent the growth rate as a function of the wavenumber. Finally, in Fig. 11, we plot the largest growth rate as a function of the temperature. The maximum value of the largest growth rate is obtained for and is given by