Phoretic interactions and oscillations in active suspensions of growing Escherichia coli

Bioluminescence imaging experiments were carried out to characterize spatio-temporal patterns of bacterial self-organization in active suspensions (cultures) of bioluminescent Escherichia coli and its mutants. An analysis of the effects of mutations shows that spatio-temporal patterns formed in standard microtitre plates are not related to the chemotaxis system of bacteria. In fact, these patterns are strongly dependent on the properties of mutants that characterize them as self-phoretic (non-flagellar) swimmers. In particular, the observed patterns are essentially dependent on the efficiency of proton translocation across membranes and the smoothness of the cell surface. These characteristics can be associated, respectively, with the surface activity and the phoretic mobility of a colloidal swimmer. An analysis of the experimental data together with mathematical modelling of pattern formation suggests the following: (1) pattern-forming processes can be described by Keller–Segel-type models of chemotaxis with logistic cell kinetics; (2) active cells can be seen as biochemical oscillators that exhibit phoretic drift and alignment; and (3) the spatio-temporal patterns in a suspension of growing E. coli form due to phoretic interactions between oscillating cells of high metabolic activity.


Introduction
Active colloids are particles that self-propel through viscous fluids by converting energy extracted from their environment into directed motion [1][2][3][4]. of artificial and natural active colloids are chemical reaction-driven Janus particles and living bacteria, respectively [1][2][3][4]. It is thought that these types of active particles differ by their mechanisms of selfpropulsion and chemotaxis. In the case of Escherichia coli and many other bacteria, the chemosensory and chemotaxis signal transduction systems enable the movement of cells towards attractants and away from repellents by controlling the rotational direction of their flagellar motors [5]. Chemotactic E. coli are capable of forming sophisticated patterns of self-organization on agar plates [6,7]. In turn, artificial Janus swimmers [8,9] are transported due to phoresis [4], which is defined as any sort of colloidal migration in gradients of all kinds [10]: solute concentration, electrical field, temperature and so on. The behaviour of self-phoretic micro-swimmers may resemble the behaviour of a living cell by exhibiting 'colloidal' chemotaxis along the self-generated chemical gradients [11][12][13][14][15][16] or gradients of fuel [17,18]. It is thought that synthetic chemotactic colloids exhibit chemotactic behaviour due to interplay between self-propulsion, phoretic drift and alignment, driven by and mediated via chemicals [12]. Moreover, the chemotaxis of synthetic colloids is sometimes modelled by chemotaxis-diffusion equations (Keller-Segel model) which are normally used to model the chemotactic behaviour of living cells [2,[11][12][13][14][15][16]19]. Recently, a representative model for phoretically interacting active colloids was introduced [19]. It was shown, in particular, that phoretic interactions induce pattern formation in typical Janus colloids when the lifetime of the phoretic field is not too small compared to the persistence time of a phoretic swimmer. Although phoretic effects have been the focus of many theoretical and experimental studies, there is a lack of insights into the role of such effects in suspensions of living cells. The complexity of transport of E. coli in liquid suspensions was revealed recently by Schwarz-Linek et al. [3].
The main goal of this paper is to understand the role of phoretic interactions in self-organization of bacteria. More specifically, we investigate whether phoretic interactions are able to enforce the formation of spatio-temporal patterns in liquid suspensions of bacteria.
It has been recently shown that non-invasive bioluminescence imaging can be applied to provide new insights into bacterial migration and self-organization in liquid cultures [20][21][22][23][24]. On the basis of the analysis of imaging data, it was hypothesized in [21,23,24] that the phoretic transport and pH-taxis of metabolically active aerobic cells may be responsible for the formation of the oscillating (or 'mergingemerging') patterns of chemotaxis (for the corresponding mathematical models, see [25][26][27][28][29]). This study is designed to elucidate the possible role of phoretic interactions in self-organization of bacteria. The paper makes the following contributions: -A number of long-run (24 h) experiments to characterize self-organization of E. coli and its mutants in standard microtitre plates are described. Mutants that potentially differ with respect to their ability to swim were chosen for the investigation. -The results of mathematical modelling of pattern formation are presented. The model is based on Keller-Segel equations of chemotaxis with logistic cell kinetics. -The underlying mechanisms of bacterial self-organization in liquid cultures of E. coli are proposed. It is argued that patterns are formed due to phoretic interactions between cells that exhibit biochemical oscillations [30,31].

Cultures and mutants
The cultures of E. coli and their mutants were prepared as described in [20,21,23,24]. Cells were grown in LB (Luria-Bertani) broth medium. To start the growth, 1 ml of overnight culture stock was added to 100 ml of LB medium. The cells were aerobically cultured in Erlenmeyer flasks on a rotary shaker at 30°C. Optical density measurements at 600 nm (OD600) were used to monitor the concentration of cells. The freshly harvested mid-log phase cultures (OD600 = 0.4-0.6) were chosen for examination [20,21,23,24]. The induction of the culture was provided by adding a stock solution of isopropyl-β-D-thiogalactoside into LB medium up to a final concentration of 1 mM. The Quick and Easy E. coli Gene Deletion Kit (Gene Bridges) was used for generation of deletion mutants. Functional cassettes for ndh, nuoA, cydA, cyoA genes flanked by homology arms were prepared by polymerase chain reaction (PCR) ( PCR products were purified by running the whole PCR sample on an agarose gel and subsequent gel extraction using the GeneJET PCR Purification Kit (Thermo Fisher Scientific). DNA concentration was adjusted to 100-400 ng µl −1 .
Transformation of E. coli DH10β (F-endA1 recA1 galE15 galK16 nupG rpsL lacX74 Φ80lacZ M15 araD139 (ara,leu)7697 mcrA (mrr-hsdRMS-mcrBC) λ-) with Red/ET expression plasmid pRedET (mediating ampicillin resistance) and disruption of a chromosomal DNA fragment by the FRT-flanked PGK-gb2-neo cassette (mediating kanamycin resistance) were performed as described in the technical protocol of Quick and Easy E. coli Gene Deletion Kit, v. 2.3. Verification of the successfully modified genome by PCR analysis was performed using oligonucleotides primer 2, primer 3, 10BndhtkrUS, 10BndhtkrDS, 10BnuoAtkrUS, 10BnuoAtkrDS, 10BtkrcydAUS, 10BtkrcydADS, 10BcyoAtkrUS and 10BcyoAtkrDS and Maxima Hot Start PCR Master Mix (Thermo Fisher Scientific), according to the manufacturer's recommendations. PCR primer2 and primer3 is designed to verify the correct insertion of the FRT-flanked kanamycin/neomycin resistance marker cassette, and the other primers are designed to verify whether the gene of interest is undisrupted, disrupted or replaced with the FRT-flanked cassette.
Other E. coli mutants were from the Keio collection [32]. All mutants were transformed by pSB417 plasmid, which harbours the luxCDABE genes from Photorhabdus luminescens placed under the control of the lac promoter [33].

Image data analysis and numerical methods
Microtitre plate-based assays were performed. Black flat bottom 96-well microtitre plates were used (Greiner Bio-one GmbH, Germany). In each microtitre plate test experiment the bioluminescence images of 12 adjacent wells (three rows, four columns) filled with 0.25 ml of liquid cultures of control and various mutant cells were acquired simultaneously using a CCD camera CoolSNAPcf (Photometrics, USA) equipped with a Schneider Kreuznach Xenon XN 0.95/25 objective. The experiments were carried out in a light-tight box at 22°C. The frames were taken once per minute for 24 h with an exposure of 1 min. The bioluminescence imaging data were collected using the image processing software Image-Pro Express V. 6.3. The freshly harvested cultures of each mutant were tested three times in 2-4 wells. Bioluminescence images (1440 images per sample) of 7-11 samples of each mutant were analysed.
Software with functions to extract pseudo-one-dimensional spatio-temporal plots from the two-dimensional images and to isolate and count high-concentration area aggregates was developed. Quantification of the self-organization process by counting the average number of unstable aggregates was performed using this software. Signal processing techniques were used to determine such an average number of aggregates. First, the bioluminescence images were pre-processed for noise reduction by a one-dimensional Gaussian filter in x-direction (each horizontal line) [34]. Second, the Otsu method [35] was used to determine a threshold for each line. With the applied threshold, the number of highlighted areas at each time step was calculated. Lazarus integrated development environment [36] and Free Pascal [37] were used to develop this software.
Owing to the nonlinearity of governing equations of the Keller-Segel-Fisher model [21][22][23][24] used, an exact analytical solution does not exist for the initial boundary value problem. Hence, the bacterial selforganization was simulated numerically using the finite-difference technique [38]. A uniform discrete grid was used for the space (x) dimension and for the time. An explicit finite-difference scheme has been built as a result of the difference approximation of the equations [38]. Verification of the scheme has been performed through varying step sizes and error tolerances. In the simulation, there were 224 grid points on the x-axis. Also a constant dimensionless step size of 0.00005 was used in the time direction. Very  [27] for simulating spatio-temporal oscillations in the Keller-Segel system. A software tool implementing computational schemes was developed using Free Pascal [37].

Patterns of bioluminescence
Patterns of inhomogeneous bioluminescence were observed in microtitre plate wells filled with cell suspensions in growth medium. Escherichia coli and its mutants show similar patterns during 24 h recordings. Illustrative snapshots are shown in figure 1a (see also electronic supplementary material, movie S1). The bioluminescence along the circular contact line and across the diameter of the well was measured using the image processing software [24]. Typical examples of the space-time plots of bioluminescence are shown in figure 1b. The patterns of bioluminescence observed during the longrun experiments can be characterized by the formation of luminous aggregates, which were also considered as azimuthal waves or plumes/clouds [20][21][22][23][24]. The number of aggregates in the microtitre plate well varies between one and six. The observed fluctuations in the number of aggregates during 24 h recordings are illustrated by figure 1c. The self-organization of bacteria in a microtitre plate well was characterized by the average number of bioluminescent aggregates (m) formed along the three-phase contact line [24].

Analysis of mutants regarding their ability to self-organize
Mutants with potentially modified chemotactic or phoretic properties were chosen for our study. These mutants lack the genes encoding various structural, regulatory or transport proteins, and enzymes, which are schematically depicted in figure 2a. The scheme was designed in accordance with the relevant literature [39][40][41][42][43][44]. The results of the analysis of mutants for their ability to form spatiotemporal patterns are summarized in figure 2b. The dynamic range of changes in the average number of aggregates (m) is rather narrow. The ratio between the largest and the smallest values of m is approximately 2. Nevertheless, a semi-quantitative comparison of the effects of mutations on the patterns of self-organization can be provided.
As can be seen in figure 2b, pattern formation is suppressed, i.e. a small average number of luminous aggregates are formed, in suspensions of mutants deficient in the proton-translocating (two protons per one electron) dehydrogenase (NDH-1) or cytochromes (CyoA, CydA) of the electron transport chain of bacteria [39]. A similar negative effect is characteristic of bacteria deficient in a major porin (OmpF) responsible for aqueous channels in the outer membrane [40]. Contrarily, mutants deficient in nonproton-translocating dehydrogenase (NDH-2) [39], flagella components (FliC or FliD) [5], the main body of fimbriae (FimA or FimI) [41], aerotaxis receptors (Aer or Tsr) [39] or proton antiporters (NhaA, NhaB or ChaA) [42] show a slightly promoted pattern formation, i.e. a higher number of luminous aggregates are formed when compared with control bacteria. The cells deficient in chemotaxis proteins (CheA, CheB or CheY) [5] flagellar motor components (MotA or MotB) [5], catalases (KatE or KatG) [43] or adhesive tips of pili (FimH or FimG) [41] generate patterns that are closest to control. Thus, the formation of bioluminescence patterns is related neither to the chemotaxis system of bacteria nor to the rotation of flagella. Seemingly, there is an alternative mechanism of the chemotactic pattern formation in suspensions of bacteria near the three-phase contact line. We explain this mechanism as 'colloidal' chemotaxis (i.e. phoretic drift and alignment of chemotactic colloids [12]), which is in line with recent studies of active suspensions of synthetic particles [11][12][13][14][15][16][17][18][19]. As follows from the analysis of the histogram shown in figure 2b, the average number of aggregates formed due to 'colloidal' chemotaxis is determined by two key factors. Firstly, it depends on the efficiency of proton translocation across membranes. The efficiency is higher (and more aggregates are formed) when the stoichiometry (H + /e − ratio) of proton translocation by the electron transport chain is higher [39] (compare the effects of dehydrogenases). The efficiency of translocation can be suppressed to some extent by hydrogen ion antiporters: mutants lacking antiporters (NhaA, NhaB, ChaA) form more aggregates than control. Secondly, the average number of chemotactic aggregates depends on the smoothness of the cell surface. 'Bald' mutants, devoid of surface appendages (fimbriae or flagella), form a higher number of chemotactic aggregates when compared with control or other mutants, which are more 'hairy'. For example, bald mutants deficient in fimbriae ( fimI) [41] tend to form more aggregates than hairy cells which lack just adhesive tips of fimbriae ( fimG or fimH). It was shown recently that surface appendages strongly impact nanomechanical and  electrokinetic properties of E. coli cells [44]. Therefore, it can be suggested that these properties may be important also in the formation of spatio-temporal patterns. In the following, we present the results of simulations of pattern-forming processes in the cultures of E. coli and discuss these processes assuming that active cells may behave as chemotactic colloids.

Mathematical modelling 3.3.1. Pattern formation in terms of the Keller-Segel-Fisher model
The Keller-Segel equations can be used to model pattern formation in bacterial populations [25]. However, similar equations can also be derived for chemotactic colloids [2,[11][12][13][14][15][16]19]. Below, our mathematical model of pattern-forming processes is provided, keeping in mind that active cells may behave as chemotactic colloids that exhibit phoretic interactions [19]. The dynamics of the E. coli population near the three-phase contact line has been described by a system of two equations, which in the dimensionless form read [21,22] and  where x and t stand for space and time, respectively, n(x, t) is the dimensionless cell density, c(x, t) is the dimensionless concentration of the chemoattractant (chemical signal), D is the diffusion coefficient, χ is the chemotactic sensitivity, α stands for the growth rate of a cell population, β stands for the saturating signal production, l is the dimensionless length of the contact line, T is the final time and is the Laplace operator. The dimensionless parameters γ and δ of the production and the degradation of the chemoattractant are usually assumed to be equal to 1 [25,27], but sometimes these parameters are kept in the dimensionless model to study the spatial pattern formation [26]. All the parameters (D, χ , α, β, γ , δ, l and T)  classification of chemotaxis models [25], these governing equations are a combination of the saturating signal production (M6) and the cell kinetics (M8) models. Governing equations (3.1) are also known as the Murray-Myerscough system [45]. The system (3.1) with zero-flux or periodic boundary conditions generates spatio-temporal oscillations [27,29,46].
At the following values of the parameters, the system (3.1) with the periodic boundary conditions produces irregular spatio-temporal oscillations qualitatively similar to the experimentally observed structures: where r is dimensionless radius of the circle corresponding to the contact line. Figure 3 shows numerically simulated space-time plots of the dimensionless cell density n(x,t), the chemoattractant concentration c(x,t) and the corresponding values n avg (t) and c avg (t) averaged on the contact line, The pattern simulation started at a 10% random initial perturbation of cells and a uniform zero concentration of chemoattractant in the medium [21,23].
As one can see in figure 3, values of n avg and c avg rather quickly (at t ≈ 10) reach the homogeneous steady-state solution (1, (γ /δ)/(1 + β)) ≈ (1, 0.58) of the system (3.1), while later n avg and c avg range notably below (0.87, 0.51), the homogeneous steady-state solution. The same effect was also noted at other values of the model parameters [22]. Additional numerical experiments assuming γ = δ = 1 showed that increasing α, D, β as well as decreasing χ leads to increasing n avg and c avg . This feature can be also revealed by a simple analysis of the mathematical model.

Pattern formation parameters
The necessary condition for instability of the homogeneous steady state was obtained by applying a linear stability analysis for equation (3.1) (see appendix A for details), This instability condition was then used to determine an interval of values for the number m of unstable aggregates, where κ = (χγ /(1 + β) 2 − Dδ − α)/2D. Note that the instability condition (3.4) is independent of the domain size, while the number m of unstable aggregates directly depends on the dimensionless radius r of the microtitre plate. Figure 4 shows the dependence of the number of unstable aggregates on the chemotactic sensitivity χ calculated at a certain cell growth rate α = 1 and on the cell growth rate α calculated at χ = 6. When calculating the numbers of aggregates, values of the other model parameters were as defined in equation (3.2).
As one can see from formulae (3.4) and (3.5), neither the instability condition nor the number of unstable aggregates depends on the cell density. As the average initial density of cells was used in defining the dimensionless parameters [22,25], the instability condition as well as the number of unstable aggregates implicitly involves also the initial density of cells. Particularly, as the dimensionless chemotactic sensitivity χ was defined as directly proportional to the cell density, all the effects of the χ -parameter can also be considered as influenced by the cell density.
The values of the dimensionless parameters of the model were determined experimentally by changing input parameters and aiming to achieve an irregular wave pattern comparable to the one shown in figure 1b. Taking into account the transformation of variables [22,25]

Effect of the chemotactic sensitivity and the cell growth
Assuming γ = δ = 1 as is usual for the dimensionless Keller-Segel-Fisher model, the chemotactic sensitivity χ and the cell growth rate α are the most important parameters contributing to the occurrence of instability of the spatially homogeneous stationary solution [26,27].

2).
A dependence (3.4) of the chemotactic sensitivity χ on the cell growth rate α is depicted in figure 5 by a solid curve and its linear approximation around a value α = 1, A similar bifurcation line in the parameter space (χ , α) was obtained when calculating hexagonal stationary solutions of a similar problem on a rectangle [27].
According to equation (3.5), the number m of aggregates varies from m min = ceil(rk 1 ) ≈ ceil(2.87) = 3 to m max = floor(rk 2 ) ≈ floor(6.86) = 6, where k 1 ≈ 1.32 and k 2 ≈ 7.53 are the unstable wavenumbers k (κ ≈ 4.52; see appendix A). The arithmetic mean of these m min and m max equals 4.5, which compares well with that (4.46) seen in the corresponding simulated pattern, marked by a solid rectangle in figure 5.

Effect of the production and degradation of the chemoattractant
To investigate the influence of the production and degradation of the chemoattractant on the formation of the spatio-temporal patterns, the effect of two dimensionless parameters γ and δ in the model (3.1) was studied.
Both the parameters, γ and δ, can be eliminated from the model (3.1) by using the following rescaling: The rescaling and dropping the asterisks yields the dimensionless governing equations the same as equations (3.1) only without the parameters γ and δ, i.e. γ and δ become equal to one in equation (3.1).
According to the rescaling (3.7), a spatio-temporal pattern simulated at any values of the parameters γ and δ can be also simulated at γ = δ = 1, while values of all other parameters are transformed using (3.7). The transformation of the chemoattractant concentration c is irrelevant when starting with a zero concentration of the chemoattractant, c(x,0), x ∈ [0,l]. An increase in the chemoattractant production rate γ can be offset only by a proportional decrease in the chemotactic sensitivity χ . The compensation of the change on the chemoattractant degradation rate δ is more complicated.
To see the relationship between the production (γ ) and degradation (δ) of the chemoattractant, the spatio-temporal patterns of the cell density were also simulated for different values of the γ -and δ-parameters at α = 1, χ = 6, keeping the other model parameters as defined in equation (3.2). The simulated patterns are shown in figure 6. The dependence derived from formula (3.4) as the production rate γ versus the degradation rate δ of the chemoattractant is depicted in figure 6 by a solid curve and its linear approximation around a value δ = 1,

Relationship between chemotactic sensitivity and chemoattractant production rate
Changing only the production rate γ and keeping all the other model parameters unchanged results in the same spatio-temporal pattern if the chemotactic sensitivity χ is also changed so that the product γ χ is kept constant.
The instability condition (3.4) can be rewritten to reveal the importance of the product γ χ for the pattern formation, According to formula (3.5), even the number of unstable aggregates directly depends on the product γ χ.
As one can see in the instability conditions given by formulae (3.4) and (3.9), the chemotactic sensitivity χ is inversely proportional to the rate γ of the chemoattractant production. For a numerical study of the relationship between χ and γ as well as the influence of the parameters on the pattern formation, the patterns were simulated by changing χ and γ while all other parameters were kept constant as defined in equation ( This expression defines the bifurcation line χγ = 5.19, which is also depicted in figure 7. As one can see in this figure, the average number of unstable aggregates near the bifurcation line varies only slightly. Four patterns marked by solid rectangles correspond to χγ = 6 > 5. 19 and are practically identical. This sameness can be simply explained by the rescaling (3.7), which becomes very simple at δ = 1. As defined in formula (3.5), the number of unstable aggregates depends on the product χγ , which is constant on the bifurcation curve χγ = const. It was already shown that the average number of aggregates at χγ = 6 is approximately equal to 4.5. This number of aggregates compares well with those seen in the simulated patterns, close to the bifurcation line shown in figure 7.
The average number of aggregates is a monotonic increasing function of χγ (see appendix A). However, figure 7 shows that the average number of aggregates slightly decreases with increasing the chemotactic sensitivity χ as well as the rate γ of the chemoattractant production. So, formula can be explained by the linearization applied to the stability analysis resulting in formula (3.5) because the linearization is effective in predicting qualitative patterns of the system behaviour only in the neighbourhood of a stationary solution.
3.4. Pattern-forming processes in the growing culture 3.4

.1. Chemotaxis and phoretic interactions
According to model (3.1), pattern-forming bacteria grow and exhibit chemotaxis. Below we provide possible physical and biological explanations of these processes in cultures of bacteria confined in the wells of microtitre plates.
The chemotactic motion of a phoretic particle can be seen as its phoretic drift and alignment due to coupling with the 'phoretic field' that is generated by all other colloids [19]. We consider a phoretic field is an effective chemical field whose gradient induces phoretic migration of cells. Thus, metabolically active bacteria, potentially, may move chemotactically (align and swim) like chemical reaction-driven Janus particles in the resulting self-produced gradient of c (phoretic field) by electrophoresis, diffusiophoresis or a similar phoretic mechanism [4,19]. The parameters χ and γ in the corresponding Keller-Segel model can be interpreted as some averages quantifying the cell surface phoretic mobility and surface activity, respectively [9,11]. Bearing in mind these remarks, it follows from formula (3.5) that the average number m of unstable aggregates formed in the microtitre plate well due to phoretic interactions can be evaluated, where r is dimensionless radius of the circle corresponding to the three-phase contact line. Equation (3.11) and the numerical values of model parameters imply that the number of unstable aggregates depends mainly on the product between phoretic mobility χ and the surface activity γ of phoretic cells. This result is consistent with the experimental data on the effects of mutations on spatiotemporal patterns of bacterial self-organization (figure 2 and the corresponding discussion). Indeed, it seems natural that 'bald' mutants, devoid of surface appendages, are more mobile than control cells or other hairy mutants, and therefore they form more aggregates. Analogously, the higher H + /e − ratio of proton translocation by the electron transport system implies a higher surface activity and higher number of unstable aggregates. Thus, the pattern-forming process in growing cultures can be associated with phoretic interactions driven and mediated by chemicals. Both the present and previous studies [21,23,24] indicate a key role of protons (hydrogen ions) in phoretic processes. It is very likely that rapidly diffusing protons are the main charged particles involved in the formation of the inhomogeneous chemical (phoretic) fields. In terms of chemotaxis, they can be considered as chemoattractants in the pH-taxis of aerobically growing and metabolically active E. coli.

Growth and oscillations
As was shown above, the dimensional cell growth rate in our model equals 0.028 s −1 . At first glance this result appears to be unexpected, because the rate of bacterial growth is of the order of h −1 [3,6,7]. In the previous studies, we interpreted similar results by assuming temporal changes of cell behaviour due to their metabolic flexibility [21,23,24]. It is known that intracellular oscillations may emerge in a living cell [30]. Despite the fact that the details of oscillatory processes, which determine the growth rate in our model, are unknown, we hypothesize that some active pole-to-pole oscillators may be involved. A well-known example of a biochemical oscillator is the Min system of E. coli [31]. This system positions the division plane of growing E. coli through pole-to-pole oscillation of Min proteins. Interestingly, the period of the swimming-resting (ON-OFF) cycle estimated above (approx. 25 s) is nearly two times shorter than the period of pole-to-pole oscillations of Min proteins in E. coli (approx. 50 s) [31]. In fact, this means that in both cases the bacterium becomes non-polar and passive every 25 s and we believe we are faced with the same oscillatory phenomenon ( figure 8). This speculative assumption should be validated in future studies.

Phoretic interactions between biochemical oscillators
In summary, the formation of spatio-temporal patterns in a growing culture can be explained by phoretic interactions in a population of biochemical oscillators. To characterize the relevant dynamic processes, let us return to equation (3.11). The terms in brackets in equation (3.11) represent the rates of processes in the pattern-forming diffusion-chemotaxis-growth system (3.1). We speculate that  Figure 8. A model of aerobically growing E. coli cell. A metabolically active cell is considered as a pole-to-pole oscillator that moves chemotactically (curved line) like a chemical reaction-driven Janus particle in the inhomogeneous chemical (phoretic) field, which is generated by all other colloids (blue triangle). The patterns of chemotaxis were observed when the lifetime of the phoretic field is much longer than the period of oscillations.
χγ /(1 + β) 2 and Dδ can be interpreted, respectively, as a growth rate and degradation rate of the inhomogeneous chemical (phoretic) field. The parameter α can be interpreted as a growth rate of the phoretic swimmer. If so, then the system can be characterized by three temporal parameters: the duration of the phoretic field formation, τ form = (1 + β) 2 /χγ , the lifetime of the phoretic field, τ phor = 1/Dδ, and the persistence time of a phoretic swimmer, τ swim = 1/α. The numerical simulations indicate that the specific oscillatory patterns form when the characteristic times (as calculated from equation (3.2)) are approximately related as follows: 0.1τ phor ∼ τ swim ∼ 2τ form . Consequently, when α = 0.028 s −1 , we have τ phor = 360 s; τ swim = 36 s; τ form = 18 s. These estimates indicate that phoretic interactions between shortperiod (approx. 0.5 min) biochemical oscillators in the growing culture of E. coli should be accompanied by the collective chemical oscillations (grow and decay of the phoretic field) with a period of about 6 min. These results show that the pattern-forming mechanisms in liquid medium and in agar plates [6,7] are different. Although, in both cases, the key processes are chemotaxis and growth, the origins and timescales of these processes are different.

Conclusion
The self-organization of E. coli cells, as detected by bioluminescence imaging, can be modelled by the Keller-Segel equations of chemotaxis with logistic cell kinetics. However, the analysis of mutants regarding their ability to self-organize shows that the observed formation of the chemotaxis patterns is related neither to chemosensory and chemotaxis signal transduction systems, nor to dynamics of flagellar motors of bacteria. This analysis together with mathematical modelling indicates that the spatiotemporal patterns in the suspensions of E. coli form due to phoretic interactions between oscillating cells of high metabolic activity.
Data accessibility. The time-lapse video is provided as the electronic supplementary information accompanying this paper.