On the stabilizing effect of chemotaxis on bacterial aggregation patterns

Based on the model by Kawasaki et al. [9], Leyva and collaborators [14] have proposed the following chemotactic, non-linear degenerate cross-diffusion system of equations in order to model the dynamics of B. subtilis on thin agar plates:

ut=∇⋅(σuv∇u)−∇⋅(ξ(u,v)χ(v)∇v)+θκuv,vt=DvΔv−uv,$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {{u_t}} \hfill & = \hfill & {\nabla \cdot (\sigma uv\nabla u) - \nabla \cdot (\xi (u,v)\chi (v)\nabla v) + \theta \kappa uv,} \hfill \\ {{v_t}} \hfill & = \hfill & {{D_v}{\rm{\Delta }}v - uv,} \hfill \\ \end{array} \end{array}$$(1)

for (x,y) ∊ Ω ⊂ ℝ², t ∊ (0, +∞). Here Ω is an open, bounded spatial domain and t > 0 denotes time. Scalar functions v = v(x,y,t) and u = u(x,y,t) represent the concentration of nutrients and the density of bacterial cells, respectively. Moreover, k > 0 is a positive constant and θ > 0 is the constant conversion factor.

The constant D_v > 0 is the diffusion coefficient for the nutrient, whereas the diffusion coefficient of bacteria is a non-linear function of the concentrations, as proposed by Kawasaki et al. [9],

Du=Du(u,v)=σuv,$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {D_u} = {D_u}(u,v) = \sigma uv, \end{array} \end{array}$$(2)

where σ > 0 measures the hardness of the agar substrate. This nonlinear diffusion coefficient was proposed to model the experimental observations by Oghiwari et al. [18], which suggest that bacteria are immotile when either the bacterial density u or the nutrient concentration v are low, and become active as u or v increase.

The chemotactic sensitivity function χ = χ(v) follows a Lapidus-Schiller receptor law [13]:

χ(v)=χ0Kd(Kd+v)2;$$\begin{array}{} \displaystyle \chi \left( v \right) = \frac{{{\chi _0}{K_d}}}{{{{\left( {{K_d} + v} \right)}^2}}}; \end{array}$$(3)

here χ₀ > 0 is a constant, measuring the strength of the chemotactic signal, and K_d > 0 is the receptor-ligand binding dissociation constant (measured in nutrient concentration units), representing the nutrient level needed for half receptor to be occupied. It depends on the bacterial species under consideration, alone. Finally, the bacterial response function ξ = ξ (u,v), which must be non-negative for attractive chemotaxis, measures the sensitivity of the bacteria to the effective sensed gradient, namely ξ(v)∇v, and it has the same units as a diffusion coefficient. In the DBM-regime of soft agar and nutrient poor medium, Leyva et al. [14] proposed the following bacterial response function:

ξ(u,v)=uDu(u,v)=σu2v.$$\begin{array}{} \displaystyle \xi \left( {u,v} \right) = u{D_u}\left( {u,v} \right) = \sigma {u^2}v. \end{array}$$(4)

This choice obeys the empiric law proposed by E. Ben-Jacob and his group [3, 8] which states that the bacterial response function should be proportional to the product of the bacterial density and its diffusion coefficient:

|ξ|∝uDu.$$\begin{array}{} \displaystyle \left| \xi \right| \propto u{D_u}. \end{array}$$(5)

The explanation by Ben-Jacob and collaborators goes as follows. Bacterial motion is, necessarily, cooperative at low bacterial densities. This can be roughly modeled by taking the diffusion coefficient to be dependent on bacterial density. Thus, the chemotactic flux should be modified accordingly as J_chem = −ξ(u,v)χ(v)∇v, where ξ follows the empirical rule (5). Notice that when the diffusion coefficient is constant, D_u(u,v) = 1 (after normalizations), then the bacterial response function is ξ = u, recovering the standard Keller-Segel [10, 11] chemotactic flux: J_chem = −uχ(v)∇v.

System (1) is further endowed with initial conditions of the form

u(x,y,0)=u0(x,y),(x,y)∈Ω,$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {u\left( {x,y,0} \right) = {u_0}\left( {x,y} \right),}&{\left( {x,y} \right)} \end{array} \in {\rm{\Omega }}, \end{array}$$(6a)

v(x,y,0)=v0,constant,$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {v(x,y,0) = {v_0},} & {{\rm{constant}},} \\ \end{array} \end{array}$$(6b)

where v₀ > 0 is the uniform initial distribution of nutrient, and the function u₀ = u₀(x,y) is the initial inoculation of bacteria in χ (for example, a localized Gaussian distribution). In addition, we impose no-flux or Neumann boundary conditions,

∇u.v=0and∇v.v=0,at ∂Ω,$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\nabla u.v = 0}&{{\rm{and}}} \end{array}}&{\nabla v.v = 0,} \end{array}}&{{\rm{at\ }}\partial {\rm{\Omega }}} \end{array}, \end{array}$$(7)

where v ∊ ℝ², |v| = 1, is the unit outer normal at each point of ∂Ω.

After an appropriate scaling of variables (see [14] for details), system (1) can be put in non-dimensional form,

ut=σ0(1+Δ^)∇⋅(uv∇v)−χ0σ0(1+Δ^)∇⋅(u2v(1+v)2∇v)+uv,vt=Δv−uv,$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {{u_t}} \hfill & = \hfill & {{\sigma _0}(1 + {\rm{\hat \Delta }})\nabla \cdot (uv\nabla v) - {\chi _0}{\sigma _0}(1 + {\rm{\hat \Delta }})\nabla \cdot (\frac{{{u^2}v}}{{{{(1 + v)}^2}}}\nabla v) + uv,} \hfill \\ {{v_t}} \hfill & = \hfill & {{\rm{\Delta }}v - uv,} \hfill \\ \end{array} \end{array}$$(8)

after suitable modifications for the functions D, ξ and χ (with a slight abuse of notation, we use the same letters for all functions and independent variables). Following Kawasaki el al. [9] (see also [14]), and in order to avoid intrinsic symmetries in the bacterial patterns induced by numerical discretization, in (8) we have considered

D(u,v)=σuv=σ0(1+Δ^)uv,$$\begin{array}{} \displaystyle D\left( {u,v} \right) = \sigma uv = {\sigma _0}\left( {1 + {\rm{\hat \Delta }}} \right)uv, \end{array}$$(9)

where σ₀ > 0 is a non-dimensional constant, and Δ^$\begin{array}{} \displaystyle \hat \Delta \end{array}$ is a stochastic fluctuation. Whence, after the scaling, the only non-dimensional free parameters are σ₀ > 0, which measures the hardness of the agar substrate (larger values of σ₀ for lower agar concentrations); χ₀ ≥ 0, which measures the strength of the chemotactic signal; and v₀ > 0, the relative initial concentration of nutrient.

We now present the results of numerical simulations of system (8), subject to boundary conditions of Neumann type (7), and with initial conditions of form (6a) and (6b). The calculations were performed in a two-dimensional square domain

Ω=(−L/2,L/2)×(−L/2,L/2),$$\begin{array}{} \displaystyle \Omega = ( - L/2,L/2) \times ( - L/2,L/2), \end{array}$$

with center at the origin, and sides of length L = 680. We employed an explicit Runge-Kutta time stepper and finite differences scheme, on a spatial square grid of N = 4096 points. Thus, the discretization mesh width was taken as Δx = Δy = L/N ≈ 0.166. The time step size considered was Δt = 0.011, for which the method is stable.

For the sake of comparison with the numerical simulations in [9, 14], the initial distribution of bacterial cell density was of Gaussian type, centered at the origin,

u0(x,y)=uMe−(x2+y2)/6.25,$$\begin{array}{} \displaystyle u_0(x,y) = u_M e^{-(x^2 + y^2)/6.25} \end{array}$$

where u_M = 0.71 is the maximum density at the center. The uniform initial distribution of nutrient was taken as

v0(x,y)≡v0=0.71,$$\begin{array}{} \displaystyle v_0(x,y) \equiv v_0 = 0.71, \end{array}$$

as considered in Kawasaki et al. [9]. The coefficient σ > 0, that appears in the diffusion for bacteria, was perturbed from a mean value σ₀ > 0 via a stochastic fluctuation Δ^$\begin{array}{} \displaystyle {\rm{\hat \Delta }} \end{array}$, following (9). The value for the non-dimensional parameter σ₀ was taken as σ₀ = 4.0, measuring the softness/hardness of the agar substrate. These parameter values (v₀ = 0.71, σ₀ = 4.0) correspond to an environment with semi-solid agar (soft medium) and low-to-medium concentration of nutrients, where DBM patterns arise, characterized by smooth envelope fronts moving outward. In the absence of chemotaxis (cf. [9,14]), these conditions produce dense branches (fingering) and ramified patterns.

While fixing these conditions for nutrient and agar substrate, we switched on the chemotactic parameter χ₀ ≥ 0. We performed simulations for different parameter values of the chemotactic sensitivity ≥₀. The results of the simulations can be observed in Figure 1. A sequence of three time steps (in rows) of the contour plots for the computed bacterial density are presented. The numerical values for v₀ = 0.71 and σ₀ = 4.0 are kept fixed in the simulations, whereas, in columns, the values for the chemotactic sensitivity are varied and taken as σ₀ = 0, 2.5 and 5.0, respectively. Notice the appearance of a colony envelope front for u, moving outward as time goes by.

Fig. 1

The first apparent change with respect to absent or weak chemotaxis is that the colony envelope front moves faster outward, as the chemotaxis is increased. This expected behavior has been already quantified in [14], where estimates on the normal speed of the front in terms of the chemotactic signal were established. The other significant change is morphological. Notice that the envelope front goes from a branched (dendrite-like) pattern (σ₀ = 0) to a more homogeneous bacterial pattern (σ₀ = 5.0), in which the protuberances seem to disappear.

Apparently, the increase in the chemotactic signal produces a transition from the DBM-morphology to the homogeneous disk type morphology which is observed in soft-agar, high-nutrient regimes [3, 9]. The presence of the chemotactic signal seems to supress the onset of instability that causes the formation of branches. This behaviour was already noticed in the simulations performed by Leyva et al. [14]. However, the authors did not provide any justification of this sudden change in morphology of the patterns.

In order to quantify the effects of chemotaxis on the branching patterns, we are going to make some approximations on the solutions to system (10). First, we employ the reduction suggested by Kawasaki et al. [9]. They observed that, due to the form of the reaction terms (production and loss terms are balanced), in the absence of diffusion the total mass is approximately conserved for small values of nutrient and bacterial density. In our model, if we neglect both diffusion and chemotaxis (the latter is also vanishing thanks to the bacterial response under consideration –hypothesis (H3)–), then we add the equations and integrate in time to obtain

u+v≈C,constant.$$\begin{array}{} \displaystyle u + v \approx C, \quad \text{constant.} \end{array}$$

This constant can be taken as C = v₀, the initial nutrient reference value. This approximation is consistent with the numerical simulations of [9] (in the absence of chemotactic signal), and those of [14], where chemotaxis is incorporated. Such simulations show that, at first order, u is constant inside the envelope front and v is near zero (total consumption of resources), whereas away from the front u = 0 and the nutrient is not yet consumed. Thus, we make the following approximation

v=v0−u.$$\begin{array}{} \displaystyle v = {v_0} - u. \end{array}$$(11)

Upon substitution into system (10) one obtains a scalar equation for the bacterial density of the form

ut=∇⋅(D(u)∇u)+g(u),$$\begin{array}{} \displaystyle {u_t} = \nabla \cdot \left( {D\left( u \right)\nabla u} \right) + g\left( u \right), \end{array}$$(12)

where the effective nonlinear, degenerate, density-dependent diffusion coefficient is given by

D(u):=D(u,v0−u)(1+uχ(v0−u)),$$\begin{array}{} \displaystyle D\left( u \right){\rm{ : = }}D\left( {u,{v_0} - u} \right)\left( {1 + u\chi \left( {{v_0} - u} \right)} \right), \end{array}$$(13)

and the effective reaction term is

g(u) := v0u(1−uv0).$$\begin{array}{} \displaystyle g\left( u \right){\rm{ : = }}{v_0}u\left( {1 - \frac{u}{{{v_0}}}} \right). \end{array}$$(14)

It is to be observed that, under the hypothesis of attractive chemotaxis towards nutrients (H3), the effective diffusion coefficient is a non-negative, density-dependent nonlinear function. Moreover, it is doubly-degenerate, in the sense that it vanishes at u = 0 and at u = v₀ (equivalently, at v = 0). In addition, notice that the resulting reaction term (14) is of Fisher-KPP (or logistic) type [6, 12].

Finally, since χ(·) is a C² function near v = 0, we approximate χ by a constant value,

χ(v)≡χ(0)=χ0≥0,$$\begin{array}{} \displaystyle \chi(v) \equiv \chi(0) = \chi_0 \geq 0, \end{array}$$

that is, we assume that the chemotactic sensitivity does not vary for v near zero. This simplification allows us to work with the parameter χ₀ ≥ 0 of the numerical simulations as a measure of the chemotactic signal, which can be switched on and off. Since for suitably normalized values of v ∊ [0,v₀] the Lapidus-Schiller chemotactic law is uniformly bounded, this approximation represents no loss of generality.

Following Funaki et al. [7], we assume that, locally, the envelope front behaves like a planar front, and that its propagation happens in one preferred direction. Thus, let us consider an approximated one-dimensional planar front solution of the form

u(x,y,t)=ϕ(x−ct)=ϕ(z),$$\begin{array}{} \displaystyle u(x,y,t) = \phi(x-ct) = \phi(z), \end{array}$$

on a strip domain,

ΩL={(x,y)∈ℝ2:−∞<x<∞,0<y<L},$$\begin{array}{} \displaystyle \Omega_L = \{ (x,y) \in \mathbb{R}^2 \, : \, -\infty < x < \infty, \; \; 0 < y < L \}, \end{array}$$

for some L > 0. Here z := x − ct denotes the (Galilean) variable of translation, and c ∊ ℝ is the constant speed of the front. Substituting the traveling front solution into (12) we obtain the following ordinary differential equation for the profile function φ,

−cϕz=D(ϕ)ϕzz+D(ϕ)zϕz+g(ϕ).$$\begin{array}{} \displaystyle - c{\varphi _z} = D\left( \varphi \right){\varphi _{zz}} + D{\left( \varphi \right)_z}{\varphi _z} + g\left( \varphi \right). \end{array}$$(15)

In pattern formation problems, it is natural to consider infinite domains and to neglect the influence of boundary conditions. Due to finite speed of propagation and the fact that we are interested in the local-in-time, local-in-space evolution near the interface of the envelope front, working on an infinite domain Ω_L means no loss of generality. Consequently, we assume that the front has asymptotic limits of the form

u±=limz→±∞ϕ(z).$$\begin{array}{} \displaystyle {u_ \pm } = \mathop {\lim }\limits_{z \to \pm \infty } \phi (z). \end{array}$$

Since the front spreads from the region occupied by the cells toward outer regions filled with nutrient, we consider

u+=0,u−=v0.$$\begin{array}{} \displaystyle {u_ + } = 0,{u_ - } = {v_0}. \end{array}$$

Moreover, we shall assume that the front is monotone, that is,

ϕz(z)<0,for all z∈R.$$\begin{equation} \phi_z(z) \lt 0, \qquad \text{for all } \, z \in \mathbb{R}. \end{equation}$$(16)

Remark 1. The existence of traveling fronts for reaction diffusion equations where the diffusion coefficient is density-dependent and degenerate, and the reaction term is of Fisher-KPP type, as well as its monotonicity properties, have been addressed in the literature by several authors; see, for example, [16, 19]. As a by-product of the existence analysis of fronts for degenerate diffusion coefficients, which decay to zero as they approach equilibria,

D(ϕ)→0,as ϕ(z)→{0,z→+∞v0,z→−∞,$$\begin{array}{} \displaystyle D(\phi ) \to 0,\;\;{\rm{as }}{\mkern 1mu} \phi (z) \to \left\{ {\begin{array}{*{20}{l}} {0,}&{z \to + \infty }\\ {{v_0},}&{z \to - \infty ,} \end{array}} \right. \end{array}$$

it turns out that the heteroclinic orbit belongs to a center manifold passing through the degenerate equilibrium points. For instance, when φ(z) → 0 as z → +∞, the center manifold on the phase plane looks like [19, 20],

ϕz=−αϕ+O(ϕ2),$$\begin{array}{} \displaystyle {\phi _z} = - \alpha \phi + O({\phi ^2}), \end{array}$$

for some constant α > 0, from which we can deduce the exponential decay of φ and its derivatives,

|ϕ|,|ϕz|,|ϕzz|≤Ce−αz,z→+∞,$$\begin{array}{} \displaystyle |\phi |,|{\phi _z}|,|{\phi _{zz}}| \le C{e^{ - \alpha z}},\quad z \to + \infty , \end{array}$$(17)

as well as the rate of decay of the diffusion coefficient

D(ϕ)=D′(0)ϕ+12D″(0)ϕ2+…=O(ϕ)≤Ce−αz,$$\begin{array}{} \displaystyle D(\phi ) = D'(0)\phi + {\textstyle{1 \over 2}}D''(0){\phi ^2} + \ldots = O(\phi ) \le C{e^{ - \alpha z}}, \end{array}$$(18)

for φ ~ 0. Similar decay rates hold at the other degenerate point, as φ → v₀ when z → -∞. See [15, 16, 19, 20] for further information.

Let us consider perturbations of the form w = w(z,y,t), satisfying the boundary conditions

w(±∞,y,t)=0,t>0,0<y<L,w(z,L,t)=w(z,0,t)=0,t>0,z∈ℝ.}$$\begin{array}{} \displaystyle \left.\begin{array}{*{20}{r}} {w\left( { \pm \infty ,y,t} \right) = 0,}&{t > 0,0 < y < L,}\\ {w\left( {z,L,t} \right) = w\left( {z,0,t} \right) = 0,}&{t > 0,z \in \mathbb{R}.} \end{array}\right\} \end{array}$$(19)

Substituting u(x,y,t) = φ(z) + w(z,y,t) into (12) and expanding the nonlinear functions we obtain

−cϕz−cwz+wt=D(ϕ)(wzz+wyy)+(D(ϕ)ϕz)z+2D(ϕ)zwz+D(ϕ)zzw+g(ϕ)+g′(ϕ)w+O(|w|2+|w||∇w|).$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} { - c{\phi _z} - c{w_z} + {w_t} = D(\phi )({w_{zz}} + {w_{yy}}) + {{(D(\phi ){\phi _z})}_z} + 2D{{(\phi )}_z}{w_z} + D{{(\phi )}_{zz}}w + }\\ {g(\phi ) + g'(\phi )w + O(|w{|^2} + \left| w \right|\left| {\nabla w} \right|).} \end{array} \end{array}$$

In view of the profile equation (15) and linearizing around the front, one arrives at the following linear equation for the perturbation

wt=D(ϕ)(wzz+wyy)+(2D(ϕ)z+c)wz+(D(ϕ)zz+g′(ϕ))w.$$\begin{array}{} \displaystyle {w_t} = D(\phi )({w_{zz}} + {w_{yy}}) + (2D{(\phi )_z} + c){w_z} + (D{(\phi )_{zz}} + g'(\phi ))w. \end{array}$$(20)

In order to establish the associated spectral problem, let us specialize (20) to solutions of the form

w(z,y,t)=eλtU(z,y),$$\begin{array}{} \displaystyle w(z,y,t) = {e^{\lambda t}}U(z,y), \end{array}$$

where λ ∊ C and U ∊ L²(Ω_L), subject to the following boundary conditions

U(±∞,y)=0,y∈[0,L],Uy(z,0)=Uy(z,L)=0,z∈ℝ.}$$\begin{array}{} \displaystyle \left.\begin{array}{*{20}{r}} {U( \pm \infty ,y) = 0,}&{y \in [0,L],}\\ {{U_y}(z,0) = {U_y}(z,L) = 0,}&{z \in \mathbb{R}.} \end{array}\right\} \end{array}$$(21)

Substituting into (20) we obtain the following spectral problem,

λU=D(ϕ)(Uzz+Uyy)+(2D(ϕ)z+c)Uz+(D(ϕ)zz+g′(ϕ))U,$$\begin{array}{} \displaystyle \lambda U = D(\phi )({U_{zz}} + {U_{yy}}) + (2D{(\phi )_z} + c){U_z} + (D{(\phi )_{zz}} + g'(\phi ))U, \end{array}$$(22)

where λ ∊ C plays the role of an eigenvalue, and U ∊ L²(Ω_L) is the associated eigenfunction. The right hand side of (22) defines a closed, densely defined operator in L²(Ω_L), namely

LU=D(ϕ)(Uzz+Uyy)+(2D(ϕ)z+c)Uz+(D(ϕ)zz+g′(ϕ))U,L:D(L)=H2(ΩL)⊂L2(ΩL)→L2(ΩL),$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {{\cal L}U} \hfill & = \hfill & {D(\phi )({U_{zz}} + {U_{yy}}) + (2D{{(\phi )}_z} + c){U_z} + (D{{(\phi )}_{zz}} + {g^\prime }(\phi ))U,} \hfill \\ {{\cal L}} \hfill & : \hfill & {{\cal D}({\cal L}) = {H^2}({\Omega _L}) \subset {L^2}({\Omega _L}) \to {L^2}({\Omega _L}),} \hfill \\ \end{array} \end{array}$$

with domain 𝒟(ℒ ) = H²(Ω_L). ℒ is the linearized operator around the front. For any U ∊ L²(Ω_L) satisfying (21) and any m = 0,1,2,..., let us define

Um(z):=∫0LU(z,y)Ym(y)dy,$$\begin{array}{} \displaystyle {U_m}(z): = \int_0^L U (z,y){Y_m}(y){\mkern 1mu} dy, \end{array}$$

where

Ym(y)={1L,m=0,2Lcos(mπyL),m=1,2,…$$\begin{array}{} \displaystyle {Y_m}(y) = \left\{ {\begin{array}{*{20}{l}} {\frac{1}{{\sqrt L }},}&{m = 0,}\\ {}&{}\\ {\sqrt {\frac{2}{L}} \cos\ (\frac{{m\pi y}}{L}),}&{m = 1,2, \ldots } \end{array}} \right. \end{array}$$

Eigenfunctions are thus decomposed as

U(z,y)=∑m=0+∞Um(z)Ym(y).$$\begin{array}{} \displaystyle U(z,y) = \sum\limits_{m = 0}^{ + \infty } {{U_m}} (z){Y_m}(y). \end{array}$$

Substituting into (22) and collecting equal modes for each m ∊ ℤ, m ≥ 0, we establish the following hierarchy of spectral equations,

λU0=D(ϕ)∂z2U0+(2D(ϕ)z+c)∂zU0+(D(ϕ)zz+g′(ϕ))U0,for m=0,$$\begin{array}{} \displaystyle \lambda {U_0} = D(\phi )\partial _z^2{U_0} + (2D{(\phi )_z} + c){\partial _z}{U_0} + (D{(\phi )_{zz}} + g'(\phi )){U_0},\quad {\rm{for }}{\mkern 1mu} m = 0, \end{array}$$(23)

and,

λUm=D(ϕ)∂z2Um+(2D(ϕ)z+c)∂zUm+(D(ϕ)zz+g′(ϕ)−m2π2L2D(ϕ))Um,for m∈ℤ,m≥1,$$\begin{array}{} \displaystyle \lambda {U_m} = D(\phi )\partial _z^2{U_m} + (2D{(\phi )_z} + c){\partial _z}{U_m} + (D{(\phi )_{zz}} + g'(\phi ) - \frac{{{m^2}{\pi ^2}}}{{{L^2}}}D(\phi )){U_m},\qquad {\rm{for }}{\mkern 1mu} m \in ,{\mkern 1mu} m \ge 1, \end{array}$$(24)

for the same eigenvalue λ ∊ ℂ. We are interested in solutions U_m ∊ L²(ℝ) to (23) and (24), for all m ≥ 0. The right hand sides of both equations are closed, densely defined differential operators of second order in L²(ℝ), with domain H²(ℝ).

Remark 2. Observe that, in view of the profile equation (15), the function Φ := φ_z ∊ H²(ℝ) is a solution to (23) with λ = 0 and m = 0. Indeed, if we differentiate (15) we obtain

0=D(ϕ)Φzz+(2D(ϕ)z+c)Φz+(D(ϕ)zz+g′(ϕ))Φ.$$\begin{array}{} \displaystyle 0 = D(\phi) \Phi_{zz} + \big( 2D(\phi)_z + c \big) \Phi_z + \big( D(\phi)_{zz} + g'(\phi) \big) \Phi. \end{array}$$

This shows that λ = 0 is an eigenvalue (associated to translation) of (22) with eigenfunction U(z,y) = Φ(z) = φ_z(z). Clearly, its expansion is given by U0=Lϕz$\begin{array}{} \displaystyle {U_0} = \sqrt L {\phi _z} \end{array}$, and

Um(z)=∫0Lϕz(z)Ym(y)dy=0,$$\begin{array}{} \displaystyle {U_m}(z) = \int_0^L {{\phi _z}} (z){Y_m}(y){\mkern 1mu} dy = 0, \end{array}$$

for all m ≥ 1, in that case.

Following [15], for each mode m ∊ ℤ, m ≥ 0, let us define the change of variables

Wm(z):=Um(z)exp(c2∫z0zdζD(ϕ(ζ))),$$\begin{array}{} \displaystyle {W_m}(z): = {U_m}(z)\exp (\frac{c}{2}\int_{{z_0}}^z {\frac{{d\zeta }}{{D(\phi (\zeta ))}}} ), \end{array}$$

where z₀ ∊ ℝ is fixed but arbitrary. Let us suppose that the eigenfunction U ∊ H²(Ω_L) ⊂ L²(Ω_L) decays sufficiently fast to zero as z → ±∞ so that, for each m ≥ 1,

Wm∈H2(ℝ).$$\begin{array}{} \displaystyle W_m \in H^2(\mathbb{R}). \end{array}$$

Remark 3. This assumption simply implies that we restrict the analysis to the set of localized perturbations that decay fast enough. At the level of the spectral problem, by a standard cut-off argument one can always approximate any generic localized perturbation by rapidly decaying ones; at the level of the evolution equations, we are simply restricting the space of controlled perturbations to a smaller space with fast decay.

Hence, we can compute

∂zUm=exp(c2∫z0zdζD(ϕ(ζ)))(∂zWm−c2D(ϕ)Wm),∂z2Um=exp(c2∫z0zdζD(ϕ(ζ)))(∂z2Wm−c2D(ϕ)∂zWm+(cD(ϕ)z2D(ϕ)2+c24D(ϕ)2)Wm),$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {{\partial _z}{U_m} = \exp (\frac{c}{2}\smallint _{{z_0}}^z\frac{{d\zeta }}{{D(\phi (\zeta ))}})\left( {{\partial _z}{W_m} - \frac{c}{{2D(\phi )}}{W_m}} \right),}\\ {\partial _z^2{U_m} = \exp (\frac{c}{2}\smallint _{{z_0}}^z\frac{{d\zeta }}{{D(\phi (\zeta ))}})\left( {\partial _z^2{W_m} - \frac{c}{{2D(\phi )}}{\partial _z}{W_m} + (\frac{{cD{{(\phi )}_z}}}{{2D{{(\phi )}^2}}} + \frac{{{c^2}}}{{4D{{(\phi )}^2}}}){W_m}} \right),} \end{array} \end{array}$$

and substitute into the spectral equation (24) to obtain

λWm=D(ϕ)∂z2Wm+2D(ϕ)z∂zWm+(H(z)−m2π2L2D(ϕ))Wm,$$\begin{array}{} \displaystyle \lambda {W_m} = D(\phi )\partial _z^2{W_m} + 2D{(\phi )_z}{\partial _z}{W_m} + (H(z) - \frac{{{m^2}}}{{{\pi ^2}{L^2}}}D(\phi )){W_m}, \end{array}$$(25)

for all m ∊ ℤ, m ≥ 0 and the eigenvalue λ ∊ ℂ associated to the perturbation. Here

H(z):=D(ϕ)zz+g′(ϕ)−c24D(ϕ)−cD(ϕ)z2D(ϕ).$$\begin{array}{} \displaystyle H(z): = D{(\phi )_{zz}} + g'(\phi ) - \frac{{{c^2}}}{{4D(\phi )}} - \frac{{cD{{(\phi )}_z}}}{{2D(\phi )}}. \end{array}$$

The following lemma states that any eigenvalue is real and non-positive, that is, the front is spectrally stable under transversal perturbations.

Remark 4. Lemma 1 means that, at first order approximation (locally planar front for a scalar equation due to the balanced source and loss kinetic terms), the envelope front is stable under transversal small, local-in-space perturbations. This behaviour is verified by the numerical calculation of the actual (curved) envelope front.

Remark 5. The result of Lemma 1 can be extrapolated to a whole family of operators indexed by the transversal Fourier frequencies. Indeed, if we consider the eigenvalue problem (22) on L²(ℝ²) and eigenfunctions of the form U = e^iξ W (z) with W ∊ L²(ℝ) then we obtain the following family of linear operators

LξW=D(ϕ)Wzz+(2D(ϕ)z+c)Wz+(D(ϕ)zz+g′(ϕ)−ξ2D(ϕ))W,Lξ:L2→L2,$$\begin{array}{} \displaystyle \begin{array}{*{20}{r,l}} {{{\cal L}_\xi }W} \hfill & = \hfill & {D(\phi ){W_{zz}} + (2D{{(\phi )}_z} + c){W_z} + (D{{(\phi )}_{zz}} + {g^\prime }(\phi ) - {\xi ^2}D(\phi ))W,} \hfill \\ {{{\cal L}_\xi }} \hfill & : \hfill & {{L^2} \to {L^2},} \hfill \\ \end{array} \end{array}$$

for each ξ ∊ ℝ. A similar analysis (under similar decay assumptions) leads to spectral stability of each operator ℒ_ξ, with ξ ∊ ℝ fixed, generalizing in this fashion the observation by Butanda [4] in the constant diffusion case.

Let us define

η(z):=D(ϕ(z))≥0,z∈ℝ,$$\begin{array}{} \displaystyle \eta (z): = \sqrt {D(\phi (z))} \ge 0,\qquad z \in \mathbb{R}, \end{array}$$

and the customary weighted energy function spaces

Hηk(ℝ;ℂ)={v:η(z)v(z)∈Hk(ℝ;ℂ)},$$\begin{array}{} \displaystyle H_\eta ^k(\mathbb{R};\mathbb{C}) = \{ v{\mkern 1mu} :{\mkern 1mu} \eta (z)v(z) \in {H^k}(\mathbb{R};\mathbb{C})\} , \end{array}$$

for k ∊ ℤ, k ≥ 0, which are Hilbert spaces endowed with the inner product (and norm),

〈u,v〉Hηk:=〈ηu,ηv〉Hk,‖v‖Hηk2:=‖ηv‖Hk2=〈v,v〉Hηk.$$\begin{array}{} \displaystyle {\langle u,v\rangle _{H_\eta ^k}}: = {\langle \eta u,\eta v\rangle _{{H^k}}},\quad \quad \left\| v \right\|_{H_\eta ^k}^2: = \left\| {\eta v} \right\|_{{H^k}}^2 = {\langle v,v\rangle _{H_\eta ^k}}. \end{array}$$

According to custom, we denote Hη0(ℝ;ℂ)=Lη2(ℝ;ℂ)$\begin{array}{} \displaystyle H_\eta ^0(\mathbb{R};\mathbb{C}) = L_\eta ^2(\mathbb{R};\mathbb{C}) \end{array}$.

Let us suppose that λ ∊ (−∞,0] is an eigenvalue associated to an eigenfunction U = ΣU_mY_m Φ H²(Ω_L) of problem (22). Let m ∊ ℝ, m ≥ 1 be a single mode for which U_m ≠ 0. Consequently, W_m ≠ 0 and ‖Wm‖Lη2>0$\begin{array}{} \displaystyle {\left\| {{W_m}} \right\|_{L_\eta ^2}} > 0 \end{array}$. (The modes m ≥ 1 for which U_m ≡ 0 do not contribute to the energy estimate.) Then from equation (29) and Lemma 1 we have that

|λ|∫−∞+∞D(ϕ)|Wm|2dz=∫−∞+∞D(ϕ)2|∂zWm|2dz+m2π2L2∫−∞+∞D(ϕ)2|Wm|2dz+J,$$\begin{aligned} |\lambda| \int_{-\infty}^{+\infty} D(\phi) |W_m|^2 \, dz &= \int_{-\infty}^{+\infty} D(\phi)^2 |\partial_z W_m|^2 \, dz + \frac{m^2}{\pi^2 L^2} \int_{-\infty}^{+\infty} D(\phi)^2 |W_m|^2 \, dz + J, \end{aligned}$$

where

J:=−∫−∞+∞D(ϕ)2Ψz∂z(|Wm|2Ψ)dz=∫−∞+∞2D(ϕ)D′(ϕ)ϕzΨzΨ|Wm|2dz+∫−∞+∞D(ϕ)2ΨzzΨ|Wm|2dz,$$\begin{array}{} \displaystyle J: = - \int_{ - \infty }^{ + \infty } D {(\phi )^2}{\Psi _z}{\partial _z}(\frac{{|{W_m}{|^2}}}{\Psi }){\mkern 1mu} dz = \int_{ - \infty }^{ + \infty } 2 D(\phi )D'(\phi ){\phi _z}\frac{{{\Psi _z}}}{\Psi }|{W_m}{|^2}{\mkern 1mu} dz{\mkern 1mu} + {\mkern 1mu} \int_{ - \infty }^{ + \infty } D {(\phi )^2}\frac{{{\Psi _{zz}}}}{\Psi }|{W_m}{|^2}{\mkern 1mu} dz, \end{array}$$

after integration by parts.

In view of last lemma we have

|J|≤C∫−∞+∞D(ϕ)|Wm|2dz=C‖Wm‖Lη22,$$\begin{array}{} \displaystyle \left| J \right| \le C\int_{ - \infty }^{ + \infty } {D(\phi )} {\left| {{W_m}} \right|^2}dz = C\left\| {{W_m}} \right\|_{L_\eta ^2}^2, \end{array}$$

for some uniform C > 0. Substituting back yields the estimate

|λ|‖Wm‖Lη22≤(supz∈ℝD(ϕ))(‖∂zWm‖Lη22+m2π2L2‖Wm‖Lη22)+C‖Wm‖Lη22.$$\begin{array}{} \displaystyle \left| \lambda \right|\left\| {{W_m}} \right\|_{L_\eta ^2}^2 \le (\mathop {\sup }\limits_{z \in \mathbb{R}} D(\phi ))(\left\| {{\partial _z}{W_m}} \right\|_{L_\eta ^2}^2 + \frac{{{m^2}}}{{{\pi ^2}{L^2}}}\left\| {{W_m}} \right\|_{L_\eta ^2}^2) + C\left\| {{W_m}} \right\|_{L_\eta ^2}^2. \end{array}$$

This estimate holds for each fixed eigenvalue λ and all values of m ≥ 1. It provides an upper bound for |λ| in terms of m (recall that |λ| = −λ > 0, in view of Lemma (1)). Notice that the effective diffusion coefficient D(φ) (that is, the weight of the Sobolev norms) depends on the intensity of the chemotactic signal χ₀ as well. Therefore, we normalize with respect to the weighted L²-norm to obtain

|λ|≤ρmCm(supz∈ℝD(ϕ))+C,$$\begin{array}{} \displaystyle |\lambda | \le {\rho _m}{C_m}(\mathop {\sup }\limits_{z \in \mathbb{R}} D(\phi )) + C, \end{array}$$(30)

where

ρm:=‖Wm‖Hη12‖Wm‖Lη22>0,$$\begin{array}{} \displaystyle {\rho _m}: = \frac{{\left\| {{W_m}} \right\|_{H_\eta ^1}^2}}{{\left\| {{W_m}} \right\|_{L_\eta ^2}^2}} > 0, \end{array}$$

and C_m > 0 is a constant depending on m, satisfying

Cm≤C1+m2C2,$$\begin{array}{} \displaystyle {C_m} \le {C_1} + {m^2}{C_2}, \end{array}$$

for some uniform constants C₁,C₂ > 0.

Therefore, estimate (30) shows that the bounds for the eigenvalues and their dependence on the intensity of the chemotactic signal are controlled by the function

supz∈ℝD(ϕ)=maxu∈[0,v0]D(u)=maxu∈[0,v0](σ0v0u(1−uv0)(1+χ0u)).$$\begin{array}{} \displaystyle \mathop {\sup }\limits_{z \in \mathbb{R}} D(\phi ) = {\max _{u \in [0,{v_0}]}}D(u) = {\max _{u \in [0,{v_0}]}}({\sigma _0}{v_0}u(1 - \frac{u}{{{v_0}}})(1 + {\chi _0}u)). \end{array}$$

Let

G(u):=u(1−uv0)(1+χ0u),u∈[0,v0].$$\begin{array}{} \displaystyle G(u): = u(1 - \frac{u}{{{v_0}}})(1 + {\chi _0}u),\qquad u \in [0,{v_0}]. \end{array}$$

After straightforward calculations one finds that, for positive values of the chemotactic sensitivity χ₀ > 0, the maximum of the cubic polynomial G on the interval [0,v₀] occurs at

u*(v0,χ0)=13(v0−1χ0+(v0−1χ0)2+3v0χ0)∈(0,v0).$$\begin{array}{} \displaystyle {u_*}({v_0},{\chi _0}) = \frac{1}{3}\left( {{v_0} - \frac{1}{{{\chi _0}}} + \sqrt {{{({v_0} - \frac{1}{{{\chi _0}}})}^2} + \frac{{3{v_0}}}{{{\chi _0}}}} {\mkern 1mu} } \right) \in (0,{v_0}). \end{array}$$

Thus, the maximum of G (and hence the energy bound for |λ| in (30)) is controlled by

Θ(v0,χ0):=G(u*(v0,χ0))=u*(v0,χ0)(1−u*(v0,χ0)v0)(1+χ0u*(v0,χ0)).$$\begin{array}{} \displaystyle \Theta ({v_0},{\chi _0}): = G({u_*}({v_0},{\chi _0})) = {u_*}({v_0},{\chi _0})(1 - \frac{{{u_*}({v_0},{\chi _0})}}{{{v_0}}})(1 + {\chi _0}{u_*}({v_0},{\chi _0})). \end{array}$$

By standard calculus tools, one easily verifies that, for fixed values of the initial nutrient concentration v₀, the bound Θ is an increasing function of the chemotactic signal χ₀ ≥ 0. Figure 2 shows this behaviour of Θ as a function of χ₀, for different fixed values of v₀ = 0.5,0.75,1.0,1.5.

Fig. 2

Whence, in view of Lemma 1 and (30), we have

0>λ>−(C+ρmCmσ0v0Θ(v0,χ0)).$$\begin{array}{} \displaystyle 0 > \lambda > - (C + {\rho _m}{C_m}{\sigma _0}{v_0}\Theta ({v_0},{\chi _0})). \end{array}$$(31)

The right hand side of (31) is a negative, decreasing function of χ₀ ∊ (0,+∞) for each fixed value of the nutrient level v₀. In lay terms, this shows that for greater values of χ₀ ≥ 0, the eigenvalues of the linearized operator around the front are allowed to be more negative. Since λ is actually the in-time growth rate of perturbations of the form w = e^λtU, this suggests that λ is further pushed to the negative side as the chemotactic signal is increased, delaying the formation of unstable patterns (fingering) around the envelope front. To sum up, although the envelope front moves faster due to chemotaxis, the spectral bounds for the eigenvalues of the linearized problem actually decrease as functions of the chemotactic sensitivity χ₀, demanding greater energy to make the patterns transversally unstable. Finally, it is to be observed that our estimation is elementary, in spite of working with a density-dependent, degenerate, nonlinear diffusion coefficient.