Travelling wave analysis in chemotaxis: case of starvation

In this paper we investigate the existence of travelling wave solutions for a chemotaxis model under the scenarios of zero growth and constant growth rate. We use Lie symmetry analysis to generate generalized travelling wave solutions, a wider class of solutions than that obtained from the standard ansatz. Unlike previous approaches, we allow for diffusivity and signal degradation. We study the influence of cell growth, diffusivity and signal degradation on the behaviour of the system. We apply realistic boundary conditions to explicitly provide biologically relevant solutions. Our results generalize known results.

. In pharmacology, chemotaxis is involved in drug delivery to the targeted defective area (Lajkó et al. 2013;Sahari et al. 2014). The beauty of the dynamics of chemotaxis is that cells manifest harmonious behaviour, while behaving independently. This was observed independently by Engelmann (1881a, b), Pfeffer (1888) and Beyerinck (1895). With the remarkable work of Adler (1966a, b) in the past fifty years, bacterial chemotaxis became one of the better-documented systems in Biology. Adler (1966aAdler ( , b, 1975 observed travelling bands of bacteria when he introduced a population of cells (E. coli) in a capillary tube accommodating oxygen and an energy source. Two bands of cells were formed; the first band consumed all the oxygen and the second band consumed the residual energy source. Bands were also observed without the adding of the energy source; cells consumed oxygen and excreted a gradient of energy source (Adler 1966b). Bak et al. (1987) noticed that the bands were in a form of a circular ring. The complexity of the geometric patterns caused by chemotaxis cannot be intuitively explained from experiments (Murray 2002). As a result, mathematical modelling approaches have been proposed which have been able to predict the geometric shape of the pattern Segel 1970, 1971a, b;Patlak 1953;Scribner et al. 1974;Hillen and Painter 2009). Segel (1970, 1971a, b) proposed, for the first time from a population-based perspective, a chemotaxis model (the K-S model) that describes the motion of slime and the formation of chemotactic bands of cells. The general form of the K-S model is written as follows: where t represents the time, b is the cell density, s the concentration of the critical substrate, χ(s) the chemotactic sensitivity, k(s) the consumption rate of substrates per cell, and µ(s) and D the diffusion coefficient of the bacteria and the substrates, respectively. Note that cell proliferation was ignored in the K-S model. A singularity in the chemotactic sensitivity was required to produce travelling wave solutions Segel 1970, 1971a;Scribner et al. 1974). Such a hypothesis is problematic, given that it can cause the bands to move with unbounded velocity (we note that the speed of the band should not be larger that the speed of a single cell) (Xue et al. 2011). This unnecessary restriction can be overcome by the consideration of other relevant factors. It was shown, for instance, that adding of logistic growth terms can lead to travelling wave solutions (with non-singular sensitivity) (Nadin et al. 2008). In the case of logarithmic chemotactic sensitivity, Wang (2013) demonstrated that the adding of substrates degradation does not produce travelling wave solutions.
From the cell-based perspective, Patlak (1953) proposed the first model for chemotaxis to depict the random walk process of a particle with external bias and persistence of direction. This model was later improved by Alt (1980) and Othmer et al. (1988). Recently, Xue et al. (2011) formulated a model which takes into account the interaction between two substrates (nutrients and attractants). What is remarkable about their (1) ∂b ∂t = ∇ · (µ(s)∇b) − ∇ · (bχ(s)∇s), (2) ∂s ∂t = D∇ 2 s − k(s)b, model is that it can be applied in a variety of biological situations, including population dynamics to describe the competition between two species from a microscopic level (the individual species behaviour). Variables describing intracellular processes such as metabolism and transduction of the signal were explicitly represented in the Xue et al. 's (2011) model. Travelling wave solutions with a unique wave speed were demonstrated in the scenario of zero growth, without requiring a singularity in the chemotactic sensitivity. We Tchepmo Govinder (2015b, 2016) extended these results by allowing for diffusivity and cell proliferation, and provided explicit solutions for the first time. Franz et al. (2014) studied the case of starvation. They assumed that cells consume chemoattractants only (which do not diffuse over the space), and considered a non constant growth of bacteria. They proved the existence of travelling wave solutions in the case of no chemotaxis. It has been proved that the parabolic limit of the microscopic model is the Keller-Segel model (Lui and Wang 2010).
In this paper, we will be looking at the individual behaviour of cells to understand the convergence and harmonization of their motion. The aggregation and movement (with constant speed) of cells are the centre of our study. We will focus on the case of low presence (or absence) of nutrients as the formation of bands of cells was observed in this situation (Adler 1966b;Brenner et al. 1998). The existence of travelling wave solutions will be investigated. Unlike previous approaches, we will allow for diffusivity, and will account for signal degradation and constant cell growth. We will also study the impact of microscale parameters (such as cell growth rate, cell unbiased turning rate and cell speed) on the macroscopic behaviour of the system. We introduce the model in "Reduced model and analysis" section. As symmetry analysis has proven to be very effective in finding useful solutions to PDEs (Clarkson 1995), we utilise that approach for our system of PDEs. We generate a class of invariants that lead to generalized travelling wave solutions. In some cases, we utilise dynamical systems analysis to further investigate the behaviour of the solutions. (This confirms our previous findings on the interplays between group theory and dynamical systems analysis Tchepmo Djomegni and Govinder 2014.) Realistic initial and boundary conditions are then applied to obtain relevant solutions. We discuss our results in "Discussion" section.

Reduced model and analysis
The model emanates from previous experiments (Blat and Eisenbach 1995;Berg 1991, 1995;Woodward et al. 1995) in which bacteria (E. coli) consume nutrients and excrete a signal gradient, then aggregate in different patterns formed in response to this gradient. We are concerned with the case of limited resources (low presence or absence of nutrients). In this scenario cells consume the excreted signal only. A set of chemical processes occur within the cells to enable them to survive and respond to their surroundings. Xue et al. (2011) developed a model to describe the intracellular metabolism, written as follows: where z = (z 1 , z 2 ) depicts the cellular metabolism, F(x, t) is the concentration of nutrients, and t f and t m are the characteristic time scales for the production of the immediate variables z 1 and z 2 , respectively. In the above description, it is assumed that after consumption of succinate F(x, t) a variable z 1 is involved to facilitate the production of the signal S(x, t) via the pathway The low level of nutrients will cause z 1 to catalytically influence the production of a starving variable z 2 , via the metabolic pathway where φ stands for the reactants/products assumed to be in excess (Xue et al. 2011). We overlook the explicit representation of the variable z 1 (given that t f ≪ t m ), and we assume fast signal transduction (given that the adaptation time of the signal transduction t a ≪ t m Xue et al. 2011). Then the distribution of the cells can be described in onedimensional space as follows: where p ± (x, z, t) is the density of cells at the position x, the internal state z = z 2 and time t, moving with constant speed ±s, and h are the turning rate function and the proliferation rate of the cells, respectively. We will consider the following turning rate function (Xue et al. 2011): where 0 is the unbiased turning rate ( 0 > 0) and χ(ξ ) = (k + |ξ |) −1 is the chemotactic sensitivity function, with k being the sensitivity coefficient. By letting and integrating (6-7) over z, one can transform (6-7) into The functions n(x, t) and j(x, t) are respectively the macroscopic cell density and the flux.
The equations describing the distribution of the signal can be given by with α and γ standing for the consumption rate and degradation rate of the aspartate, respectively. Note that the case k = 0 corresponds to unbounded sensitivity to the signal (Xue et al. 2011). As a result, the function 1 becomes the switch function The situation of no chemotaxis (i.e., χ = 0) takes place when k → ∞, and we have 1 = 0 . Here, cells are not sensitive to the signal. Our analysis will focus on these two limiting cases of unbounded sensitivity to the signal and no chemotaxis. Thought it is very hard to obtain travelling wave solutions technically if k varies, it has numerically been shown that increasing chemotactic responses results to an increase in the wave propagation (Franz et al. 2014).
The system (11-13) was analysed in the case of no chemotaxis by Franz et al. (2014). They assumed that D S = γ = 0, and h(S) is a linear function of S. Due to the complexity of the system, diffusivity has always been ignored in the mathematical analysis. It is important to note that diffusivity plays a stabilizing role in the behaviour of the system (Rosen 1977). As a result, in our analysis, we will allow for diffusivity, and will investigate the existence of travelling wave solutions under zero growth and constant growth scenarios [We note that demonstrating the existence of traveling wave solutions is equivalent to demonstrating the existence of solutions to (11-13) (Lui and Wang 2010)]. The impact of the growth rate in the behaviour of the solutions will be explored.

Lie symmetry analysis
Examining the interplay between group theory and stability analysis, we found (Tchepmo Govinder 2016, 2014) that the Lie symmetry analysis can generate new types of solutions (unlike the standard travelling wave ansatz) that play a significant role in the stability of the system. A partial differential equation of order n, as a symmetry if (Bluman and Anco 2002) where ξ i (x, y) and η j (x, y) are the infinitesimals of the Lie group of invariant transformations of (15), and G [n] is the nth extension of G (Bluman and Anco 2002). The operator G (15) E(x, y, ∂y, . . . , ∂ n y) = 0, in (16) helps to reduce the number of independent variables and the order of equations. We note that a linear combination of symmetries to an equation is also a symmetry to that equation. The expression of G [n] can be obtained in Bluman and Anco (2002).
In the case of zero growth, using (19-20), the characteristic equations associated with G are (refer to Bluman and Anco 2002) These lead to the new invariants We note from (27), that travelling wave solutions can exist, with c being the speed of the wave. As we require that solutions should not blow up as x → ±∞ or t → ∞, we take c 3 = c 4 = 0. Therefore, from the definition of the flux j (see (9)), we obtain c 1 = 0.
In the case of constant growth with α 0 � = 2 0 , the characteristic equations associated with G, now using (21-22), are and lead to the invariants As before, we take c 1 = c 3 = c 4 = 0 for physically viable (bounded) solutions. The same conditions apply in the case of α 0 = 2 0 . As a result in our analysis we use the following invariants: in all cases. We note that (37-38) are generalized travelling wave solutions (Polyanin and Zaitsev 2004;Tchepmo Djomegni and Govinder 2016). The coefficients c 2 can produce to damped (when c 2 < 0) or growing (when c 2 > 0) solutions. The case c 2 = 0 lead to the standard ansatz travelling wave solutions. Then the system (11-13) can be rewritten in term of the new invariants as follows: where the superscript ′ stands for the total derivative with respect to u. When c = s, (40-42) can be reduced to a system of two equations in three unknowns. Choosing J(u) and S 1 (u) to depend on N(u), we can demonstrate travelling wave solutions for a constant distribution of N(u) [by simply solving the second order ODE with constant coefficients (42)]. The analysis for the case of Poisson distribution (or normal distribution via asymptotic analysis) of N(u) is similar to the analysis when c � = s. We will focus in the rest of this work on the case of c � = s (with 0 < c < s). Inspired by the numerical investigations of Xue et al. (2011), we will be looking for solutions admitting a single peak of S. We note that this restriction is less important when k → ∞. Hence we assume S 1 ∈ Y S , where In our context, travelling wave solutions n(x, t) and S(x, t) must be positive, continuous and bounded, with S 1 ∈ Y S .

Case of zero growth
Here h(S) = 0. Assuming that N(u) and J(u) decay to zero as u → ∞, (40-42) can be reduced to In the case of high chemotactic sensitivity (the limiting case k → 0), the turning rate function becomes a switch function and corresponds to unbounded sensitivity to the signal (Xue et al. 2011). For S 1 ∈ Y S , the solution n(x, t) is given by where σ 1 = 2 0 /(s + c) and σ 2 = 2 0 /(s − c). Here, the total cell population is given by T = sN (0)/ 0 (obtained by integrating N(u) over the whole line R).
We assume D S = 0, then S 1 (u) is positive and continuous. When γ + c 2 ≥ 0, S 1 (u) is monotonically increasing (S ′ 1 (u) > 0). As a result, S 1 ∈ Y S cannot hold. Since S ′ 1 (u) is continuous, for S 1 to hold in Y S when γ + c 2 < 0, it is necessary that zero must be the only extremum point of S 1 (u) (the maximum), with αN (0) = −(γ + c 2 ) (because S ′ 1 (0) = 0). Then, for u < 0, Again S 1 / ∈ Y S . Non-diffusing travelling wave solutions with a single peak of S do not exist.
In the case of diffusivity, substituting (47) into (46) and integrating, we obtain , the coefficients c i j are constants of integration, and the functions I k i (v) and K k i (v) are the two linearly independent solutions to the modified Bessel's equation.
Proof We choose c 2 ≤ 0 in order to produce damped solutions. We note that the functions I k i (v) and K k i (v) given in (49)  This implies that and converges to zero as u → −∞ (provided γ + c 2 > 0). Then S(x, t) converges to zero as x → ±∞ or t → ∞ (with x − ct < 0). Now we assume u ≥ 0. As u → ∞, which implies that as u → ∞. Therefore, if c 2 2 � = 0 , S(x, t) will blow up as x → ∞ (with t fixed). However, if c 2 2 = 0, S 1 (u) converges to zero as u → ∞ (given that I k 2 α 2 e −(σ 2 /2)u converges to zero u → ∞). As a result, S(x, t) converges to zero as x → ∞ or t → ∞.
Given that the only point in which S 1 (u) may not be continuous is zero, and that S 1 (0) exists, then the convergence of S(x, t) at the boundaries guarantees the boundedness of S(x, t).
Proof Invoking Proposition 1, we note that n(x, t) and S(x, t) are positive, continuous and bounded, and the function e −(c/(2D S ))u is monotonically decreasing. We only need to show that S 1 ∈ Y S . When u ≥ 0, the function I k 2 α 2 e −(σ 2 /2)u is decreasing. Therefore, S 1 (u) is also monotonically decreasing, and converges to zero as u → ∞.
When u < 0, the function I k 1 α 1 e (σ 1 /2)u is increasing, and from (50), we have Given that the function f (u)e −(c/(2D S )u is monotonically increasing, then S 1 (u) is also monotonically increasing. Consequently, S 1 ∈ Y S .
In the case of no chemotaxis (i.e., χ = 0), we have 1 = 0. Then from (45), N(u) blows up as u → −∞; cells can only aggregate in the half plane u ≥ 0. Travelling wave solutions satisfying S 1 ∈ Y S do not exist. This is consistent with Xue et al. (2011). However, if we relax the assumption on S 1 , we can demonstrate travelling wave solutions.
Theorem 2 In the absence of chemotaxis, generalized travelling wave solutions for the system (11-13) exist (without the restriction S ∈ Y S ) and are explicitly given in the case of non diffusivity (D S = 0) by and S(x, t) = S 1 (u)e c 2 t , where with u ≥ 0, σ 0 = 2 0 c/(s 2 − c 2 ) and c 2 ≤ −γ, and in the case of diffusivity (D S � = 0) by □ The proof follows the same procedure as Theorem 1. The solutions are illustrated in Fig. 2.
Assuming α 0 � = 2 0 , N(u) and J(u) are given by where For α 0 < 2 0 , 1 and 2 have the same sign (we note that 1 2 = −α 0 (α 0 − 2 0 )/(s 2 − c 2 ) ). Then S 1 / ∈ Y S , because bounded solutions will be represented only in a half plane. When α 0 > 2 0 , then 1 > 0 and 2 < 0. We will choose C 1 and C 2 so that N(u) will not blow up as u → ±∞. For u < 0, we will take C 2 = 0, and for u ≥ 0, we will take C 1 = 0 . This will require discontinuity of the flux at zero. Non-diffusing travelling wave solutions (admitting a single peak of S) do not exist. However, the requirement that S 1 ∈ Y S is less important in the case of no chemotaxis. We remark that Xue et al. 's (2011) results reflected this relaxation. In our case we do find travelling wave solutions (see Fig. 3).
Theorem 3 Non-diffusing (D S = 0) generalized travelling wave solutions for the system (11-13) exist without the restriction S 1 holding in Y S (with χ = 0). They are explicitly given by and S(x, t) = S 1 (u)e c 2 t , where for α 0 = 2 0 (with −γ ≤ c 2 < 0 and u ≤ 0). In the case of α 0 � = 2 0 , the solutions are given by γ 4 e 2 u , u ≥ 0, and S(x, t) = S 1 (u)e c 2 t , where for α 0 > 2 0 (with c 2 = −γ), and by (63) and S(x, t) = S 1 (u)e c 2 t , where for α 0 < 2 0 (with u ≤ 0 and −γ ≤ c 2 < 0, if 1 > 0, or u ≥ 0 and c 2 = −γ, if 1 < 0). □ The non-diffusing solutions S 1 (u) of Theorem 3 are obtained directly by integrating the first order system (59-61). We note that all of the solutions are bounded, for they are continuous and converge at the boundaries. A negative flux (see (66)) simply means that most of the cells move to the left (recall that j = s(n + − n − )). Unlike Franz et al. 's (2014) results, we do not require a minimal wave speed. This therefore constitutes a generalization of their findings.
We note that the discontinuity of the flux at zero does not necessarily imply N(u) to be discontinuous at zero. In fact, for J (u) = sN (u) and the initial conditions given by Theorem 4 For D S � = 0, α 0 > 2 0 , −γ < c 2 ≤ 0 and γ 3 J (0 − ) = γ 4 J (0 + ) = N (0), travelling wave solutions (with S 1 ∈ Y S ) for the system (11-13) exist and are explicitly given by The proof of the above theorem is similar to the proof of Theorem 1. See Fig. 4 for an illustration of the solutions.

High chemotactic sensitivity (k → 0)
Here 1 is given by (14). For S 1 ∈ Y S , the system (40-42) becomes where for u > 0, and γ 4 e 2 u , u ≥ 0, for u < 0. When α 0 = 2 0 , the coefficient b 1 vanishes and we obtain We note that the total cell population is given by T 2 = 2sN (0)/(2 0 − α 0 ), and the flux can be negative. This simply means that most of the cells move to the left [we recall that j = s(n + − n − )]. The analysis in this situation is mathematically similar to the case of zero growth in Sect. 2.2. Travelling wave solutions satisfying S 1 ∈ Y S exist only in the case of diffusivity, underlying the importance of introducing this biological process into the model.
When at most one of the constants C i is zero (for instance C 2 ), Eq. (77) is difficult to solve explicitly for S 1 (u) (given the form of N(u) when u > 0). In this situation, we only look at the asymptotic behaviour of the solutions as u → ±∞. We note that the origin is the only equilibrium point of the system (75-77), and the determinant of the corresponding Jacobian matrix around the origin is given by We consider α 0 < 2 0 to guarantee the stability of N(u) and J(u)(see (83)). The eigenvalues dictating the behaviour of S 1 (u) are given by For c 2 + γ < 0, 11 and 22 are both negative; non growing solutions are possible only in the half plane u ≥ 0. Then S 1 cannot hold in Y S . For c 2 + γ > 0, 11 and 22 have opposite signs. To obtain convergence, we will choose the initial data so that the behaviour of S 1 (u) will be controlled only by 22 and 11 as u approaches +∞ and −∞, respectively (this method was also applied in §2.3.1). In this situation S 1 (u) is positive, since none of the eigenvalues is complex (We note that the nonlinear term does not affect the stability of S 1 , given that the eigenvalues are nonzero.). The challenge in getting explicit solutions for S 1 (u) prevents us from checking whether the restriction S 1 ∈ Y S holds for all real u γ 0 e 2 u , u < 0, Djomegni SpringerPlus (2016) 5:917 or not. However, this restriction can be guaranteed as u → ±∞. In fact, from (82), as u → ±∞, where µ 1 = min( 1 , 2 ), µ 2 = min( 3 , 4 ), and δ 1 , δ 2 are coefficients of the dominant terms e µ 1 u and e −µ 2 u . Substituting (90) into (77) then integrating, and taking into consideration the boundedness conditions, one obtains, as u → ±∞, , and c 1 1 and c 2 1 are positive constants. We previously proved that S 1 ∈ Y S . Asymptotic travelling wave solutions are possible. The asymptotic behaviour of S 1 (u) is illustrated in Fig. 7.

Discussion
In this paper we studied the existence of travelling wave solutions (with a single peak for the signal) of a microscopic model for chemotaxis. We focused on the case of starvation; cells in this situation consume signal only. The effect of microscale parameters in the stability of the system was examined. Unlike previous approaches, we allowed for degradation of signal (γ � = 0). While we will compare our results to those previously obtained, it must be borne in mind that this important biological process was not considered in other results. We performed a Lie symmetry analysis to generate a large class of invariants leading to generalized travelling wave solutions. Only relevant invariants were considered, but we believe that rich information could have been extracted from the full form of invariants in different contexts. We provided explicit solutions, many for the first time.
We first considered the case of zero growth. Such a scenario is possible if the time interval is shorter than the period required for cell proliferation. When we imposed no chemotaxis, we could not find travelling wave solutions satisfying S 1 ∈ Y S . We note that Xue et al. (2011) also indicated the absence of travelling wave solutions when there is no chemotaxis. However, their results held in the case of no diffusivity (D S = 0). We have shown here that these results also hold in the diffusing case. However, when we relaxed the restriction on S 1 (less important in the absence of chemotaxis), we obtained both diffusing and non-diffusing travelling wave solutions, distributed in a half plane (see Fig. 2).
If we now consider the high chemotactic limit, we find that non-diffusing (D S = 0) travelling wave solutions do not exist. This is in contrast to Xue et al. 's (2011) results in which they were found to exist. The degradation of the signal removed this possibility in our results. However, as evidenced in Theorem 1, the incorporation of diffusion does allow for the existence of travelling wave solutions. It is interesting to note that diffusion, in a sense, counteracts the wave eradication effect of the degradation of the signal.
In order to model the behaviour of cells over more realistic time frames, we incorporated cell growth into our model. As a first attempt we assumed constant growth. In general, with no chemotaxis, non-diffusing (D S = 0) travelling wave solutions do not exist. However, we observe that this occurs due to the requirement that S 1 ∈ Y S . This is not necessary in the case of no chemotaxis. Relaxing this restriction leads to non-diffusing travelling wave solutions (see Theorem 3). Note that, unlike Xue et al. (2011), we dot not require a minimal wave speed. We are also able to find diffusing travelling wave solutions (with a discontinuous flux) provided the growth rate dominates the dynamics (see Theorem 4).
If we consider the high chemotactic limit we find that non-diffusing travelling wave solutions do not exist. However, incorporating diffusivity leads to the possibility of travelling wave solutions (see Theorems 5 and 6 ). Note that these are the first results in the case of high chemotaxis with non-zero growth. In contrast to Segel's (1970, 1971a) results in the macroscopic model under zero growth, none of our travelling wave solutions required a singularity in the chemotactic sensitivity.
When cells are highly sensitive to signals, allowing for diffusivity, we observe for α 0 < 2 0 that the total cell population (given by T 2 = 2sN (0)/(2 0 − α 0 )) increases as the growth rate α 0 becomes large; here most of the new born cells remain in the band. For α 0 = 2 0 , the growth rate controls the behaviour of the system (see § 2.3.2). In this case, the total cell population (given by T 3 = 2sN (0)/α 0 ) decreases as α 0 is large. This is due to the local depletion of the signal; some cells will move towards regions with higher concentrations of signal (this is typical in chemotactic systems). The aggregated cells here do not disperse as we demonstrated the existence of travelling wave solutions in this situation. However, for α 0 > 2 0 , we notice instability. The cell growth rate controls the behaviour of the system, and prevents formation of the aggregation.
The inverse phenomenon is observed in the limiting case where cells are not sensitive to the signal gradient; they move randomly in this situation. For α 0 ≤ 2 0 , we obtained instability in the system, travelling wave solutions do not exist. However, for α 0 > 2 0 , the stability of the solutions is controlled by the growth rate α 0 . We imposed restrictions on the initial conditions in order to foster a collective behaviour. Travelling wave solutions then resulted. This result is in agreement with Lauffenburger et al. 's (1984) findings (in the macroscopic model), in which travelling wave solutions exist due to the balance of growth, death and random motility.
As result of our investigation, we remark that cell growth and cell unbiased turning rate play an important role in the stability of the system and the aggregation of cells. We also remark that the total cell population in the case of zero growth (T 1 ) is less than that of the case of constant growth (T 2 and T 3 ). The wider band of cells is obtained in the case of no chemotaxis (we recall that here the growth rate α 0 controls the stability of the system, the absence of sensitivity to stimuli keep most of cells in the band). The distribution of cells are displayed in Figs. 1a, 2a, 3a, d, 4a, 5a and 6a. We also note that the total cell population T i decreases as 0 becomes large; the permanent change of direction does not necessarily destabilize the formation of bands of bacteria.

Conclusion
We have shown that it is crucial to consider the individual response of cells when studying their macroscopic behaviour. This helps us to capture microscale information which play a significant role in the system. Our results can be summarized in Table 1. For future work, we will consider a higher dimensional space and will investigate the geometric shape of bands of bacteria.