Travelling waves in hybrid chemotaxis models
Hybrid models of chemotaxis combine agent-based models of cells with partial differential equation models of extracellular chemical signals. In this paper, travelling wave properties of hybrid models of bacterial chemotaxis are investigated. Bacteria are modelled using an agent-based (individual-based) approach with internal dynamics describing signal transduction. In addition to the chemotactic behaviour of the bacteria, the individual-based model also includes cell proliferation and death. Cells consume the extracellular nutrient field (chemoattractant) which is modelled using a partial differential equation. Mesoscopic and macroscopic equations representing the behaviour of the hybrid model are derived and the existence of travelling wave solutions for these models is established. It is shown that cell proliferation is necessary for the existence of non-transient (stationary) travelling waves in hybrid models. Additionally, a numerical comparison between the wave speeds of the continuum models and the hybrid models shows good agreement in the case of weak chemotaxis and qualitative agreement for the strong chemotaxis case. In the case of slow cell adaptation, we detect oscillating behaviour of the wave, which cannot be explained by mean-field approximations.
Keywords:hybrid model travelling wave bacterial chemotaxis
The wavelike spread of cell populations plays a fundamental role in many biological processes, including development Landman:2007:MEI (), wound healing Witte:1997:GPW () and tumour invasion Gerisch:MMC:2010 (). Bacterial populations show similar phenomena, with the pioneering studies of Adler Adler:1966:CB () confirming the capacity of an E. coli population to form travelling bands via chemotaxis to extracellular signals. Mathematically, the extent to which chemotaxis can generate and sustain stationary travelling bands has motivated a number of studies, including the Keller-Segel model of Adler’s experiments which is written in the form of coupled partial differential equations (PDEs) Keller:1971:TBC (). This early model necessitated a biologically unrealistic singularity in the chemotactic sensitivity to generate stationary travelling waves: a requirement that allows bacteria behind the wave to acquire infinite speeds and to avoid “dropping-out”, an effect that leads to gradual dispersal of the band Xue:2011:TWH (); Franz:2012:HMI ().
This singularity requirement can be circumvented by incorporating other processes. The well known Fisher’s equation Fisher:1937:WAA () demonstrates travelling waves in systems coupling diffusion with logistic growth terms Fisher:1937:WAA (). Parabolic chemotaxis models with non-singular sensitivities but incorporating either logistic Landman:2003:CCM (); Landman:2005:DCD (); Nadin:2008:TWK () or non-logistic Kennedy:1980:TWS (); Satnoianu:2001:TWN () growth terms also admit travelling wave solutions. Other studies have shown that introduction of more complex nutrient terms can give rise to travelling waves, even when growth is absent Saragosti:2010:MDB (); Saragosti:2011:DPC (). An experimental system which also included two chemicals – a chemoattractant and a nutrient source – was presented in Budrene:1991:CPF (); Budrene:1995:DFS (), with stationary or transient travelling waves obtained according to t he formulation of the model Brenner:1998:PMC (); Xue:2011:TWH (). Travelling waves in chemotaxis models have also been recently studied in Li:2012:SPW (); Li:2011:ANS (); we also note the articles Horstmann:2004:UPK () and Wang:2013:MTW () for a review and analysis of travelling waves in PDE-based models. A comparison between mesoscopic (hyperbolic) and macroscopic (parabolic) PDEs has been presented in Lui:2010:TWS ().
Relatively little exploration has been conducted into travelling wave formation for chemotactic models extending beyond PDE systems, in particular those introducing terms to account for the inherent noise of biological systems. One exception is the study of Chavanis:2010:SKS (), in which a multiplicative noise term was introduced into the Keller-Segel model and the existence of travelling waves has been demonstrated within this setting. Hybrid models, in which an individual-based model for bacterial behaviour is coupled to a continuum description of extracellular signals, naturally introduce stochastic effects and will be the focus of the present paper. Such a hybrid model was formulated in Franz:2012:HMI () where it was shown that under finite cell speeds only transient travelling waves formed, even with singular chemotactic sensitivity. The individual-based model was formulated in terms of the velocity-jump model with internal dynamics Erban:2004:ICB (); Erban:2005:STS (); Xue:2009:MMT () and, in this paper, we extend the model in Franz:2012:HMI () to incorporate proliferation and death of bacteria. We analyse this system numerically and analytically with respect to its travelling wave properties, employing the biologically inspired chemotactic sensitivity presented in Xue:2011:TWH () and a linear growth term. We show that stationary travelling waves can be observed even in the absence of chemotaxis, although wave speeds are substantially increased in its presence.
The organisation of the paper is as follows: the full hybrid model is presented in Section 2 along with illustrative simulation results, while the corresponding continuum equations are derived under certain assumptions in Section 3; in Section 4 these continuum equations are analysed with respect to travelling wave properties; in Section 5 where a computational analysis and comparison of the models is presented; finally, we discuss our observations in Section 6.
2 Hybrid model of bacterial chemotaxis
In this section we formulate the hybrid model of bacterial chemotaxis which will be investigated in this paper. The model is motivated by the behaviour of the bacterium E.coli and, in its most general form, includes cell movement, sensing and response to a chemical signal, consumption of the chemoattractant, cell proliferation and death. However, for analytical tractability, we will also explore simplified hybrid models which exclude some of these processes. Bacteria are modelled as agents with internal dynamics that represent the signal processing and response of each individual while the extracellular chemical is modelled using a PDE to describe its spatio-temporal concentration. The mathematical framework and simulation techniques are reviewed in Franz:2012:HMI (). We consider the model in an effectively one-dimensional domain representing a long but narrow tube, similar to the experimental set up considered in Adler:1966:CB ().
The motion of E. coli bacteria is controlled through the coordinated rotation of flagella distributed over the cell surface Berg:1975:HBS (). Counterclockwise rotation generates a propulsive bundle that results in straight line motion of the bacterium – a so-called “run” Berg:1972:CEC (). Alternatively, clockwise rotation results in the outward flaying of flagella and a “tumble” – rotation with insignificant displacement. At the end of each tumble the bacterium chooses a new direction of movement, seemingly at random, and returns to the run phase. The lengths of the individual phases are independent from each other and distributed exponentially, yet they can be influenced by internal dynamics Berg:1975:HBS ().
Internal dynamics of the E. coli bacteria possess two principal features Bourret:1991:STP (): a quick excitation phase followed by slower adaptation. Specifically, changes in the extracellular signal concentration lead to quick excitation of the internal metabolism, signified through altered chemical concentrations inside the cell. Following excitation the internal concentrations revert slowly to normal in an adaptation process, even when the external signal remains at the raised level.
2.1 Velocity jump model with internal dynamics
Run-and-tumble dynamics are aptly modelled as a velocity-jump process Othmer:1988:MDB (); Erban:2004:ICB (). We denote by the number of bacteria (agents) in the system at time . The current state of the -th individual, will be described using its position , its velocity and a set of internal state variables that represent the states of components in the intracellular signal transduction network.
Here we concentrate on a cartoon version of the internal dynamics of bacteria written in terms of two internal variables Othmer:1998:OCS (); Erban:2004:ICB (), i.e . Internal variables and are governed by the equations
where is the excitation time, is the adaptation time, and is the concentration of chemoattractant at the position of the bacterium at time . Furthermore, bacteria move with the velocity governed through a velocity jump process with a turning frequency that depends on the internal dynamics. In this paper, we will use the biologically motivated nonlinear turning kernel developed in Xue:2011:TWH (). Hence, the full model of one individual over (a small) time step can be written as:
where and are positive constants.
In addition to the behaviour of an individual bacterium we define a signal-dependent proliferation function . We thereby interpret a positive value of as a proliferation rate, meaning that in the infinitesimal interval a bacterium at position generates an exact copy of itself with probability . Similarly, a negative value of means that the bacterium disappears (dies) with the probability . In this paper, we will use the following form for the proliferation rate :
where and are positive constants.
2.2 Evolution of the extracellular chemoattractant
For the extracellular signal we formulate a PDE that incorporates diffusion (with diffusion constant ) and signal consumption by bacteria, the latter with signal dependent rate . The equation for therefore takes the form
For the remainder of the paper we employ a linear form for the consumption function :
where is a positive constant.
2.3 Illustrative example
The hybrid model framework presented in Sections 2.1 and 2.2 includes essential features of the more complicated hybrid chemotaxis models formulated in Erban:2005:ICB (); Xue:2011:RSS (). In this section we numerically show that these processes can give rise to travelling waves. For the numerical simulation we employ techniques described in Franz:2012:HMI (). In particular, for the extracellular signal , this means that the simulation is performed on the one-dimensional domain with initial condition and zero-flux boundary conditions. We consider regularly spaced grid points , where and the values of are advanced by a small time step and a forward Euler update rule:
In the above is the symmetric, normalised and non-negative kernel
where the kernel width is a positive real number. Here, represents the influence a bacterium at position has on grid point .
The simulation of the individual bacterium is given in the full system (2.2)–(2.6) and complemented by the birth and death processes described in Section 2.1, where we use the same time step as in (2.10). To calculate the necessary off-grid values of extracellular signal, we linearly interpolate from the two nearest grid points. We further simplify the system (2.2)–(2.6) by exploiting the separate time scales for excitation and adaptation (i.e. ): specifically, we assume the update equation (2.5) for is in a quasi-equilibrium, which is identical to the assumption . The value for can therefore be calculated by
Illustrative results are presented in Figure 1. For this simulation, bacteria were initialised at positions , randomly generated as the absolute value of a Gaussian random variable with variance much smaller than the domain length . The initial velocity (direction of movement) is generated uniformly at random and initial values of the extracellular signal and internal variables are taken as
where . We simulate the system until time and plot both the distribution of bacteria and concentration of chemoattractant in Figure 1(a). We also estimate the wave speed as a function of time in Figure 1(b).
We clearly see formation of a travelling band of bacteria, moving rightwards with average speed (plotted as the dashed line in Figure 1(b)).
Influence of the growth term
To investigate the influence of the growth term on the existence of travelling waves, we simulate the full hybrid model (2.2)–(2.6) and (2.10) including () and excluding () growth and death processes. We use identical parameters to those described above and present the results in Figure 2. In Figure 2(a) the position of the wave front (defined as the right-most position for which ) is compared. The full hybrid system (dashed line) generates a straight line, indicating a wave moving with constant speed. While the system excluding growth and death (solid line) moves with a similar initial speed, speed is gradually lost over time: the shape of at different times for this case is shown in Figure 2(b). We clearly see that no true travelling wave forms, with many agents being left far behind the wave front, leading to its slowing down. Thus, we can interpret growth and death terms in terms of a stabilising role on the wave profile: although not all agents can keep up with the wave, new agents are constantly created at the front and the agents that drop out eventually die, resulting in a travelling band of agents.
3 From hybrid models to macroscopic PDEs
In this section we derive macroscopic PDEs for the spatio-temporal density of bacteria at given position and time . An implicit assumption of the derivation is spatial independence of bacteria, which allows formulation of a continuous mesoscopic system. We then use results from Erban:2004:ICB () to obtain the macroscopic equations. To illustrate the successive formulation of models we construct two systems of PDEs – denoted System (A) and System (B) – to be referred to in the remainder of the paper.
3.1 System (A)
We define the mesoscopic densities for left and right-moving bacteria, depending on their position , their internal variable and . If the signal profile was uninfluenced by bacteria, densities would satisfy the following system of hyperbolic PDEs:
The signal dynamics is described by (2.8) which can be rewritten in terms of as
The system (3.1) (for the one-particle distribution) can be derived by integrating the probability distribution function for the many particle system, utilizing the fact that the movement of individuals are biased by the signal function , but independent to each other. However, for the hybrid chemotaxis models described in Sections 2.1 and 2.2, individual bacteria interact via the extracellular signal which complicates the derivation of (3.1). In Erban:2012:ICB (), a kinetic description has been derived for a model of interacting locusts, using a modified version of the BBGKY hierarchy from the classical kinetic theory of gases Cercignani:1994:MTD (). The system we consider here is much more complicated to analyse than the locust model studied in Erban:2012:ICB (), due to the variable number of bacteria and internal variables. Thus the kinetic description (3.1) can only be considered as an approximation to the one particle distributions of the interacting system.
The capacity of the above mesoscopic system to generate travelling bands analogous to those observed in the hybrid model is illustrated in Figure 1(c)-(d). For details of the numerical method employed for this and other simulations of the continuous model, we refer to Xue:2011:TWH (). The qualitatively and quantitatively close correspondence in solutions under equivalent parameters and initial conditions corroborates the use of the above approximation.
3.2 System (B)
We consider a macroscopic model in this section. Define the macroscopic densities
and let them satisfy the following system
where the turning rates are given by
We will denote (3.5) and (3.7) along with the definition of in (3.6) as System (B). According to the analysis in Erban:2004:ICB (); Xue:2009:MMT (), System (B) is quantitatively consistent with System (A) when the external signal changes slow enough such that cells are close to their fully adapted state, in which case cell movement is only moderately modified by the signal.
In the rest of the paper, we assume diffusion of extracellular signal to occur on a much slower time scale than the active motion of the bacteria, hence . The number of parameters of the above models can be reduced by setting to one through rescaling. We show this in detail for System (B) as follows. Rescaling the variables , , , and the parameters , , taking (2.7) and substituting into System (B) we obtain, after dropping hats for notational simplicity,
We are interested in travelling wave solutions that develop from a pointwise inoculation of cells into a domain containing uniformly distributed nutrient . In this scenario, (defined as in each system) should form travelling pulses while forms a travelling front and relevant boundary conditions will be
Note that is currently unknown; we determine its value in the travelling wave analysis of Section 4. Since and are physical quantities, we search for nonnegative travelling wave solutions, i.e.
It is clear that a travelling wave of this form cannot exist for (extinction of bacteria) or for (infinite growth) and we will therefore only consider systems that satisfy . In the next section we analyse System (B) with respect to travelling wave solutions in order to obtain further insight. To do that, we use the rescaled system (3.8).
4 Travelling wave analysis
In this section we first apply the standard travelling wave ansatz to system 3.8 and derive a necessary condition for the existence of non-negative travelling wave solutions. We then reduce the resulting ODE system to two components through a change of variables and utilizing an invariant manifold identified for the problem. Finally we use phase plane methods to analyse the existence and properties of travelling wave solutions.
4.1 A necessary condition for the existence of travelling wave solutions
where the primes denote derivatives with respect to the travelling wave variable . Note that any point on the -axis is a steady state of the system 4.1 and that linear stability of such a steady state, , is governed by the eigenvalues of the matrix , where
The eigenvalues of are
Under the boundary conditions 3.9 we look for nonnegative solutions to 4.1 connecting steady states and . To admit such a solution the latter must be a stable node, since a stable spiral would imply negative values for . Hence, a necessary condition is , which is equivalent to
Given it is easy to show that .
A necessary condition for the existence of nonnegative travelling wave solutions of the system is
The above condition is reasonable, as we expect the run duration to occur on a much faster time scale than proliferation processes.
4.2 Dimension reduction
Let us now perform a change of variables by introducing the cell density and the cell flux . The travelling wave system 4.1 can then be written as
where the boundary conditions for this system are
Integrating and applying the boundary conditions at , an invariant manifold of the problem is given by
For , we obtain
where we chose the solution to the quadratic equation for that satisfies the boundary conditions as .
It can be easily shown that has two solutions in the region for all as follows. Since , is monotonically decreasing for and monotonically increasing for . With , this implies and, using for , we obtain the existence and uniqueness of the second root of : we call it . The existence of and the negativity of for , together with the condition , implies that as given in 4.11 is positive everywhere, and that the given solution therefore satisfies the nonnegativity condition.
4.3 Steady states and their linear stability
Using the two roots of and under the condition (4.5), it is clear that there are two steady states of the system 4.9-4.10: and . Linearising the system 4.9-4.10 about its steady states generates a system of the form
where, for the general steady state , we have
The eigenvalues of are identical to as given in 4.2. The steady state is therefore a stable node for all with as defined in 4.4. Similarly, it can be seen that the steady steady is a saddle point. The eigenvectors corresponding to the eigenvalues take the form
In the plane, the slopes of the eigenvectors are given by
For the steady state this slope can be written in the form
where we define similarly to (4.3).
4.4 Case I: No chemotaxis ()
We first consider the case where the chemotactic sensitivity (given by (3.6)) vanishes, i.e cells do not respond chemotactically to changes in . Here, travelling waves are generated solely through proliferation of bacteria at the wave front. To understand the wave behaviour we perform a phase plane analysis for the ODE system (4.9)–(4.10). Using (i.e. ), it reduces to
Thus, the slope of a trajectory in the plane can be written as
Additionally, an expression for the nullcline is given by
and the -nullcline is simply
Let us now show that travelling waves exist for the reduced system 4.15.
For the case (which is equivalent to ), a unique travelling wave solution for the system exists for all .
is an invariant region of the system (3.8). Since is non-decreasing everywhere in and is non-negative for and , we need only to show that the direction field on the segment points from the top half of the plane above this segment towards the bottom. Since is strictly increasing we require
where we used (4.14) in the first step and the relation for all . Using the fact that and are negative, we can use the definition of and the fact that to obtain
where we used throughout the derivation. We can therefore conclude that is an invariant region of the system (3.8). Noting that at the steady state the unstable manifold has a positive slope (), i.e. it points into the region , and using the fact that is strictly increasing inside for we can conclude that, for each , there is a heteroclinic orbit starting from and finishing at , corresponding to a travelling wave solution of the PDE system (3.8).
4.5 Case II: Increasing chemotaxis ()
It is noted that the above slope is larger than that for the non-chemotaxis case within the region of interest . Due to this increase the region for the proof of Theorem 1 is no longer invariant for this system and a travelling wave solution to (3.8) does not necessarily exist for all . The -nullcline for the full ODE system (4.9)–(4.10) is given as the solution of the quadratic equation
For a given wave speed , the -nullcline can therefore be calculated as
We can see that as due to its leading order term . Therefore, as becomes large, no -nullcline exists and is positive everywhere. Additionally, might have further roots and, in particular, might be negative in parts (or the whole) of region . This again means that is strictly growing in these parts of the domain.
We detect three different types of behaviours of trajectories starting close to , plotted in Figure 4. In particular, we can see each of these behavioural types for different values of and despite different configurations of the nullclines. In the top two plots of Figure 4 we present the case of a diverging solution. Examining ODE (4.9), we observe that for large , grows quicker than and the divergence can be identified as a finite-time blow-up. In the second case, depicted in the two plots in the middle of Figure 4, the trajectory converges to the steady state , but does so after entering the region and thereafter the region . Note that the steady state is still a stable node in this case and that this overshoot is therefore not a spiralling effect. Since these trajectories do not correspond to a non-negative solution of the ODE system (4.9)–(4.10), they do not represent travelling wave solutions to the original problem. The last case, presented in the plots on the bottom of Figure 4, corresponds to an acceptable solution and is characterised by the convergence to without crossing the line .
4.6 Case III: Infinite chemotactic sensitivity ()
As decreases further we observe that the minimal wave speed necessary to allow a non-negative travelling wave solution of (3.8) increases. In the limit , the ODE system (4.9)–(4.10) no longer has convergent solutions. However, in this limit the linearisation assumption leading to these ODEs and the system (3.8) is no longer valid and we must consider the original turning kernel as defined in (3.2). In the limit the turning rate in the hybrid model therefore becomes
Hence, bacteria moving in a favourable direction do not turn, indicating that the wave speed achieved in this limit should evolve to . In Xue:2011:TWH () it was shown, for a slightly different turning kernel, that travelling waves can exist even without growth terms and that their wave speed satisfies .
5 Computational analysis of the wave speed
In this section we computationally compare wave speeds from the hybrid model with those of the fully continuous models. Specifically, we investigate the regimes in which the latter provide an acceptable insight into the travelling wave behaviour of the hybrid model, and where they differ. We begin by investigating the non-chemotaxis case, where the minimum wave speed for the continuum systems was determined in (4.4). In Section 5.2 we show how the wave speed depends on the value of , and correspondingly the chemotactic sensitivity in the macroscopic model. A comparison with hybrid models without cell proliferation is given in Section 2.3. We conclude this section with a discussion into the effect and origin of oscillations observed under increasing the adaptation time .
5.1 Case I: No chemotaxis ()
In Section 4.4 we analysed the macroscopic PDEs in the absence of chemotaxis. Travelling wave solutions were shown to exist for all wave speeds , with determined by (4.4). In Figure 5(a), variation of (4.4) as a function of is illustrated; we note that wave speeds determined through simulation of the PDE systems correspond exactly (to accuracy of the numerical approximation) with the analytical wave speeds. We now numerically investigate the wave speed for the case in the hybrid model.
For our simulations we consider the same parameters and methods as described in Section 2.3: specifically, we set the system parameters , and . For the computations we consider a time step , a spatial resolution of on a domain with length , and simulate the system until the value of at falls below . The profiles at this time, together with the time when at falls below , are used to estimate the wave speed.
The measured wave speed for varying is illustrated in Figure 5(a), along with as predicted from the travelling wave analysis. While the relationship is similar in shape, we note that at all values of tested the measured wave speed lies below the analytical value . In the literature it has been observed that inaccuracies in numerical schemes can lead to an increase in wave speeds Reitz:1981:SNM (), therefore rendering the lower wave speed seen in Figure 5(a) as counter intuitive.
Nevertheless, we can provide the following explanation for the distinct values in the continuum and hybrid models. For the zero-chemotaxis case, wave generation and movement is solely determined by growth ahead and death behind the wave. In the continuum model an outermost “fractional bacteria population” can extend significantly beyond the wave front, since some proportion of the initial population never turns left, and hence far into the region where is very close to its initial value of . Yet this fractional population still grows exponentially (), seeding the growth and expansion of the population. The finite/discrete nature of the hybrid model precludes any fractional bacterium: the forward “tail” is necessarily finite and growth will not occur beyond the outermost individual.
For the above explanation to hold we would expect a dependence of the measured wave speed on the initial number of bacteria : continuous densities provide a closer approximation under larger numbers of bacteria and we would expect convergence in the wave speed to . Simulations in Figure 5(b) demonstrate this property, corroborating our interpretation.
5.2 Case II: Increasing chemotaxis ()
In the second set of numerical experiments we measure the dependency of the wave speed on the critical parameter , i.e. we determine the effect of increasing chemotaxis as decreases. We compare the results measured for the hybrid system with the continuous Systems (A) and (B).
We use the same parameters as in Section 5.1 and results are shown in Figure 6. The results demonstrate the regimes where correspondence across the varying modelling levels occurs: while the hybrid model (dotted line) corresponds well with its closest continuous version (mesoscopic System (A), red solid line) over a wide range of , it only corresponds with System (B) (black dashed line) for larger , diverging as decreases. Note that the turning rate (3.6) used for System (B) becomes negative at small values of and we limit the range of studied accordingly.
At larger all three models converge to a value close to as grows: in this regime the main assumption proposed for the linearisation () holds and we obtain good quantitative agreement. While this assumption becomes less acceptable as we decrease , leading to the divergent behaviour described above, we note that all models show the same qualitative agreement: increasing chemotactic responses leads to an increase in the wave speed. Note that the results for System (B) can be identically replicated using the ODE system (4.9)–(4.10) and a search algorithm for the smallest value of that admits a nonnegative solution to the system.
These numerical experiments demonstrate that chemotaxis has a significant effect on the speed of movement and that the waves cannot solely be explained by growth and death terms. Rather, we interpret birth and death processes as stabilisers to what would otherwise be transient waves Franz:2012:HMI (); Xue:2011:TWH (). This interpretation is in agreement with the results presented in Figure 2, as the initial wave speed for the system without growth seems to be similar to the wave speed of the system including growth and death terms.
5.3 Oscillations in the wave speed
An additional observation we made during the numerical experiments of the hybrid model is that