Jeans type instability for a chemotactic model of cellular aggregation

Jeans type instability for a chemotactic model of cellular aggregation

P.H. Chavanis Laboratoire de Physique Théorique, Université Paul Sabatier, 118 route de Narbonne 31062 Toulouse, France
11email: chavanis@irsamc.ups-tlse.fr
To be included later
Abstract

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

1 Introduction

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 ():

(1)
(2)

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 ():

(3)
(4)
(5)

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.

In order to make the connection between the parabolic model (1)-(2) and the hyperbolic model (3)-(5), we consider a model of the form

(6)
(7)
(8)

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 ():

(9)

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

We consider an infinite and homogeneous stationary solution of Eqs. (6)-(8) with , and such that

(10)

Linearizing the equations around this stationary solution, we get

(11)
(12)
(13)

where we have set . Eliminating the velocity between Eqs. (11) and (12), we obtain

(14)

Looking for solutions of the form and , we get

(15)
(16)

where . These equations have non-trivial solutions only if the determinant of the system is equal to zero yielding the dispersion relation

(17)

The condition of marginal stability () corresponds to and the condition of instability is

(18)

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

When the chemical has a large diffusivity, the temporal term in Eq. (8) can be neglected jager (). Thus, we formally consider . In that case, the dispersion relation (2.1) reduces to

(19)

The discriminant is

(20)

and the two roots are

(21)

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

(22)

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

(23)

If this condition is fulfilled, the unstable wavenumbers are determined by

(24)

The growth rate of the perturbation with wavenumber is . Therefore, the most unstable mode is the one which maximizes . It is given by

(25)

The largest growth rate is then determined by

In particular, for , we have

(27)

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 ).

Figure 1: Graphical construction determining the range of unstable wavenumbers. We have taken , and (solid line) or (dashed line).
Figure 2: Evolution of the growth rate of the perturbation as a function of the wavenumber for different values of . We have taken , and .

2.2.2 If :

In that case, . The unstable wavenumbers are determined by

(28)

The most unstable mode is and the largest growth rate is given by

(29)

2.2.3 If :

In that case, the system is unstable for

(30)

The growth rate diverges when

(31)

corresponding to . Close to the critical wavenumber , we have

(32)

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

(33)

Taking the limit , we find that

(34)

which is finite for but diverges like when . For and , the dispersion relation can be simplified in

(35)

For we recover Eq. (32) and for we recover Eq. (34). For , we can easily express as a function of according to

(36)

On the other hand, for , Eq. (2.1) reduces to

(37)

The positive root of this equation is

(38)

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).

Figure 3: Graphical construction determining the range of unstable wavenumbers. We have taken , and .
Figure 4: Evolution of the growth rate of the perturbation as a function of the wavenumber. We have taken , and .

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

(39)

and the two roots are

(40)

where

(41)

Writing the solution in the form

(42)

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

(43)

and stable otherwise. The range of unstable wavelengths is determined graphically in Fig. 5.

Figure 5: Graphical construction determining the range of unstable wavenumbers. We have taken , and (solid line) and (dashed line).

2.3.1 If :

In that case, the system is unstable for

(44)

and stable otherwise. This is the same criterion as for the inertial model. A necessary condition of instability is that

(45)

If this condition is fulfilled, the unstable wavenumbers are such that

(46)

These results are unchanged with respect to the inertial case.

Figure 6: Evolution of the growth rate of the perturbation as a function of the wavenumber for the Keller-Segel model. We have taken , and .

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

(47)

with

(48)
(49)
(50)

There will be two roots and provided that and . This last condition can be written

(51)

with

(52)
(53)
(54)

The discriminant of Eq. (51) is given by

(55)

The condition to have two roots and is equivalent to with

(56)

(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.

(57)

2.3.2 If :

In that case the system is unstable for the wavenumbers

(58)

The growth rate of the perturbation as a function of the wavenumber is plotted in Fig. 7. As discussed in Sec. 2.2.3, the case is special and will be considered specifically in the next section.

Figure 7: Evolution of the growth rate of the perturbation as a function of the wavenumber for the Keller-Segel model. We have taken , and .

2.4 The case and

If we neglect the temporal term in the Keller-Segel model (), we obtain the dispersion relation

(59)

so that is explicitly given by

(60)

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

(61)

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

(62)

This divergence is regularized if . Taking , i.e. in Eq. (39), we get

(63)

For , we obtain

(64)

which is finite for but diverges like when . For and , Eq. (39) can be simplified in

(65)

For we recover Eq. (62) and for we recover Eq. (64). The solution of Eq. (65) is

On the other hand, for , Eq. (39) leads to

(67)

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

(68)

On the other hand, for , taking (i.e. ) in Eq. (2.1), we get

(69)

Finally, for and , we have

(70)

which reproduces the correct behaviours (68) and (69).

3 Analogy with the Jeans problem in astrophysics

3.1 The damped Euler equations

Let us consider a particular case of Eqs. (6)-(8) corresponding to and where and are arbitrary functions. In that case, the hydrodynamical equations take the form

(71)
(72)
(73)

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

(74)

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

(75)
(76)
(77)

When , these equations are isomorphic to the damped Euler-Poisson system describing self-gravitating Brownian particles gt (); virial (). In the strong friction limit, we get

(78)

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

We consider an infinite and homogeneous stationary solution of Eqs. (75)-(77) with , and such that

(79)

For a pure Newtonian interaction with , we have . Linearizing the equations around this stationary solution, we get

(80)
(81)
(82)

where we have introduced the velocity of sound . Eliminating the velocity between Eqs. (80) and (81), we obtain

(83)

Looking for solutions of the form , we get

(84)
(85)

From these equations, we obtain the dispersion relation

(86)

In the case and , we recover the usual Jeans dispersion relation jeans ():

(87)

3.3 Instability criterion

If we set , the dispersion relation becomes

(88)

The two roots are

(89)

with

(90)

Accordingly, the system is unstable if

(91)

and stable otherwise. A necessary condition of instability is

(92)

If this condition is fulfilled the unstable wavelengths are such that

(93)

The wavelength which has the largest growth rate is given by

(94)

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

(98)

For , one has

(99)

3.5 Isothermal gas

Figure 8: Graphical construction determining the range of unstable wavenumbers.
Figure 9: Maximum wavenumber and most unstable wavenumber as a function of the temperature . The line determines the separation between stable and unstable states.

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

(100)

the growth rate of the perturbation can be written

The condition of instability reads

(102)

and a necessary condition of instability is . For the unstable wavenumbers (see Fig. 8) are such that with

(103)

The wavenumber with the largest growth rate is given by

(104)

and the largest growth rate by

(105)

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

(106)
Figure 10: Growth rate of the perturbation as a function of the wavenumber. We have taken and .
Figure 11: Dependence of the largest growth rate with the temperature. We have taken .

For