Bioconvective Dynamics: Dependence on Organism Behaviour

Bioconvection, one of the earliest reported collective behaviours of micro-organisms (Wager, 1911), refers to the spatial patterns that develop in suspensions of various swimming microbes (see Pedley and Kessler, 1992), which tend to propel themselves upwards and have a mass density greater than that of water. Upward swimming is usually a response to an external force field such as gravity or a biochemical stimulus such as a gradient in oxygen concentration. The process of pattern formation begins when the accumulation of ‘heavy’ suspended organisms near the top of the fluid generates an unstable density profile in the suspension. The excess density of the suspension is the product of the volume fraction of organisms and the difference between their density and that of water. The gravitational instability that ensues creates structures characterized by scales much greater than the size of the organisms. Patterns develop from descending blobs or curtains containing high cell concentrations, turning into plumes, then convection rolls and finally convection cells.

The stage for bioconvection is set by the breaking of up/down symmetry as a result of the differentially oriented swimming behaviour of micro-organisms. Unlike collective behaviours that originate in chemical, electrochemical or tactile signalling of cells or micro-organisms to each other, the directional swimming behaviour at the root of the bioconvection phenomenon can derive from external influences that tend to induce oriented swimming of each individual, regardless of what others do. Examples of such influences are gravity and illumination. An internally generated influence is the molecular gradient that results from consumption and unidirectional supply of specific molecular species. Although this field is collectively generated, it does not involve communication in the usual sense. The actual long-term dynamics of the circulation is determined by fluid mechanics and by the upward swimming of the organisms, required to supply the gravitational potential energy needed to compensate for viscous dissipation.

One might infer that biology and physics are separable in the mechanism of bioconvection: swimming behaviours of individual organisms, coupled only indirectly by consumption-generated gradients or gravity, is the biology; hydrodynamics, acting on the density inhomogeneities created by swimming is the physics. However, this inference is only partly true. Once the bioconvective dynamics activates itself, it alters the biological environment as well as the physical one, through mixing, long-range transport, shading and the generation of hydrodynamic flows that orient and convect. Bioconvection is therefore an exciting, complex, yet experimentally tractable, approach for studying the mutual interdependence of physics and biology, where the overall phenomenology greatly exceeds its primitive components (Kessler, 1989; Kessler and Hill, 1997). The distinguishability of microbial species by their bioconvection patterns is a very significant result that lends motivation to the experiments and to the analysis described here.

Classes of bioconvective systems can be distinguished on the basis of the mechanism of directional motion, or taxis, of the cells. (i) Some micro-organisms swim upwards because they are bottom-heavy (geotaxis or gravitaxis). (ii) The oxygen concentration gradient, generated by the oxygen consumption of the cells and by supply from the air interface, can induce upward swimming towards regions of higher oxygen concentration (aerotaxis, oxygen taxis or, in general, chemotaxis). During fully developed bioconvection, oxygen-charged water is convected downwards from the air/water interface, so the bacteria swim towards a complex, convection-dependent, usually not vertical, gradient. (iii) In phototaxis, organisms swim towards or away from a source of light. Self-shading of concentrated organism populations is a type of consumption (of photons), and is analogous to the consumption of oxygen molecules in aerotaxis.

Childress et al. (1975) proposed the first extensive theory for bioconvection of gravitactic swimmers. Their equations included both directional and stochastic aspects of the swimming velocity, but did not take into account the effect of local vorticity on swimming direction. A continuum theory for bioconvection of suspensions of aerotactic bacteria that swim upwards in molecular gradients created by consumption must include three sets of coupled equations: (i) the continuity and Navier–Stokes equations, driven by variations in the concentration of the micro-organisms that correspond to locally averaged variations in the ‘external’ driving force density, (ii) the equation for conservation of organisms, including directional and stochastic aspects of swimming together with advection by the motion of the water in which they are suspended, and (iii) the equation for conservation of the molecules that the organisms consume, including the rate of consumption by the cells, diffusive supply and advection. This set of coupled nonlinear equations and the methods for modelling the part of the organism flux that is due to swimming have been discussed by Kessler and Wojciechowsi (1997) and Kessler et al. (2000). A full solution of these equations has not been attempted. Various approximate approaches have been developed for all the above-mentioned types of bioconvection (gravi-, gyro-, photo-and aerotactic): through linear and weakly nonlinear stability theory (Harashima et al., 1988; Pedley and Kessler, 1992; Vincent and Hill, 1996; Hillesdon and Pedley, 1996; Metcalfe and Pedley, 1998; Bees and Hill, 1999; Ghorai and Hill, 2000) and through closed-form solutions of severely idealized versions (Kessler, 1985, 1986). Unfortunately, there has not been much quantitative analysis of experiments available for testing theoretical developments. Two recently published systematic quantitative studies concerned the gravitactic patterns of the algae Chlamydomonas nivalis (Bees and Hill, 1997) and the onset time of (oxygen) chemotaxis-driven bioconvection in cultures of Bacillus subtilis (Jánosi et al., 1998).

In this paper, we analyze the geometry and time-dependence of bioconvection patterns due to aerotaxis of B. subtilis. We employ methods similar to those used by Bees and Hill (1997) for C. nivalis, and we compare our results with theirs. This set of controlled experiments is the first to explore quantitatively how the bioconvection patterns of chemotactic bacteria depend on the average cell number density (ρ) in the culture and the depth (h) of the suspension. These two control parameters are the most relevant ‘physical’ ones available under fixed environmental conditions. Fourier analysis and autocorrelation were used to determine the dominant wavelength λ₀ and spectrum at the onset of the instability, for various values of ρ and h. We also measured the increase in λ₀ with time during the evolution of fully developed bioconvection. Surprisingly, we found that, although the onset of bioconvection resulting from either gravitaxis or chemotaxis relates generically to the Rayleigh–Taylor instability (Plesset and Winet, 1974), quantitative measurements revealed that both the dependence of the patterns on the control parameters and the evolution over time of the scale of the convection cells are completely different in the two cases.

The experiments were performed with the motile strain YB886, a derivative of the naturally transformable Bacillus subtilis strain 168 (Yasbin et al., 1980). This bacterium is a gram-positive, endospore-forming, aerobic eubacterium commonly found in soil.

Liquid cultures of YB886 were grown in GM1/GM2, Spizizen’s minimal glucose salts medium (Spizizen, 1958) supplemented with 0.1 % yeast extract, 0.02 % casein hydrolysate, 5 mmol l⁻¹ MgCl₂ (GM1) or 0.5 mmol l⁻¹ CaCl₂ (GM2), and the appropriate amino acids at a final concentration of 50 µg ml⁻¹. Bacteria were proliferated in GM1 at 37 °C with continuous shaking (250 revs min⁻¹). The cell number density was monitored using a Klett–Summerson colorimeter with filter no. 66.

Before the experiments, bacteria were kept in GM1 for 90 min after reaching the end of the exponential growth phase, as described by Kessler and Wojciechowski (1997). The suspensions were then diluted tenfold in warm GM2 and incubated at 37 °C with aeration for a minimum of 60 min.

Experiments were performed by transferring concentrated suspensions to Petri dishes (diameter 8–12 cm) up to a depth h. The initial cell concentration was 10⁻⁸ to 10⁻⁹ cells cm⁻³ (see below and Table 1). The system was kept at room temperature (26 °C), and no nutrients were added to the medium (GM2) during that time. A homogeneously mixed initial configuration of dissolved oxygen and bacteria was prepared by gently swirling the Petri dish for 10 s. The typical pattern formation time was of the order of minutes. The swirling/waiting cycles were repeated several times in a given dish without replacing the suspension.

The cell division rate is slow (30 min or more) compared with the formation of patterns (under 3 min). The cell concentration can therefore be considered to be constant during any particular pattern-formation process. However, the cell concentration does not remain constant over the typical length of an experimental run (1–2 h). To follow the changes in the cell population, the optical density (OD) of the suspension was repeatedly measured by turbidometry (Bausch & Lomb Spec 21). A sample of the mixed suspension was transferred to a cuvette, and its optical density was measured at 600 nm. The statistical error involved with the optical density values was estimated from repeated read-outs of the samples after shaking. Standard cell-counting procedures revealed that an optical density of 0.85 units corresponds to 6×10⁻⁸±2×10⁻⁸ cells cm⁻³ and that the relationship between optical density and cell density is linear.

Video recordings were taken of the Petri dish with dark-field illumination, as described by Jánosi et al. (1998). The recordings were then digitized using a Silicon Graphics workstation with a built-in video-capture board. The resolution of the digitized pictures was 640 pixels×480 pixels with 256 grey levels. The NTSC system used provides 30 images s⁻¹, which determined the maximum sampling rate.

The starting time of an experimental cycle was defined by observing the cessation of visible macroscopic motion left over from mixing. The average intensity before the development of optical nonuniformity, during the initial 30 images of each run, i.e. during the first second, was used to create a reference background. A sampling rate of 1 frame s⁻¹ was used for subsequent image processing, which is permissible because the characteristic time for pattern formation is much longer than 1 s. Difference images were created by subtracting the grey level values of the pixels of the background obtained at the beginning of each run. This process eliminated any problems arising from inhomogeneities in the illumination.

Except for the data illustrated in Fig. 3, pattern analysis was performed at the time of maximum contrast of the images. Maximum contrast is a convenient and reproducible marker for the onset of the instability (Jánosi et al., 1998). Briefly, the time of maximum contrast can be determined by following the evolution over time of the difference between the maximum and minimum greyness levels of the digitized images. At the beginning, small seeds of sinking plumes appear as small bright spots on a grey background. As they grow, the channels between the plumes become more and more depleted of cells, resulting in increasing contrast. A later decrease in the contrast occurs when the plumes hit the bottom of the container and the cells spread back into the depleted channels. Fourier and correlation techniques were used to extract the dominant wavelengths of the bioconvection patterns.

The Fourier method was implemented as in Bees and Hill (1997) with the following steps. (i) A central region of 256 pixels×256 pixels was cut out from each image to fit into the discrete Fast Fourier Transform (FFT) algorithm (Press et al., 1992). A typical image is shown in Fig. 1A. (ii) The Hann window (Press et al., 1992) was used to eliminate the effect of sharp edges. This windowing weights the information in the centre of the image (see Fig. 1C) and removes oscillatory errors from the spectrum. (iii) Two-dimensional Fourier transforms were performed using the FFT method. This procedure yields a complex array containing the amplitude and phase information. Fig. 1B shows the resulting amplitude pattern for Fig. 1A (note that the phase information is also necessary for the inverse transformation shown in Fig. 1C). (iv) The Fourier intensity spectrum I(k) was calculated for each pattern as:

The inset in Fig. 1D shows G(r) calculated for the pattern of Fig. 1A. The first (negative) minimum gives the characteristic distance between the brightest and darkest regions, and the first positive maximum gives the dominant wavelength λ₀.

With the appropriate calibration, both λ₀ (in mm) and the wavenumber k (in mm⁻¹) were measured. The autocorrelation function and the power spectrum mutually determine each other for sufficiently large and homogeneous data (Wiener–Khinchin theorem). In our case, the size of the samples is limited by the video resolution, and the assumption of spatial homogeneity is not strictly satisfied. Thus, both methods can be considered as approximate, and the difference between the characteristic lengths obtained by them can be used as a good measure of their inaccuracy. In Fig. 2, a correlation plot is shown for the dominant wavenumber λ₀. The scatter of the values gives the overall wavelength estimation error Δλ=±0.15 mm. The resolution of the wavelength determination from the pair-correlation function is 1 pixel (=0.162 mm). The apparently finer resolution of the FFT method is a consequence of the simple shape of the spectra, which allows a good fit to the centre of the peak.

The experiments reported and analyzed in this paper concern the evolution over time of bioconvection of the oxygen-consuming, oxygen-tactic bacteria Bacillus subtilis. Recordings traced the development of patterns, starting from initially well-mixed bacterial cultures. Fig. 3 illustrates a typical time evolution from the beginning of the instability, through the development of sharp ‘curtains’ that break up laterally, into plumes that spread at the floor of the Petri dish (Fig. 3b), to hexagons, the final fully developed state (Fig. 3g). The dominant wavelength λ₀ increases approximately logarithmically over time, as demonstrated in Fig. 3h.

Dominant wavelengths were extracted from the digitized video recordings by Fourier (i) and correlation (ii) methods. The cell concentration was determined by measuring the optical density of the sample concurrently with the recordings. Data characterizing the onset of the convection in 49 experiments, performed with various control parameter values, are summarized in Table 1.

In Bees and Hill (1997), the function I(k) was fitted by an unnormalized double Gaussian distribution. Remarkably, we found that the spectra on a double-logarithmic plot of I(k) versus k have a pronounced tent shape (Fig. 4), indicating power-law decays for both small and large wavenumbers. On the basis of this observation, the fitting function used to extract the leading wavenumber was:

We found that the cell concentration (optical density) was the most important factor determining λ₀ (Fig. 5). In fact, the empirical data are best fitted with a linear regression including both control parameters λ₀:

Bacillus subtilis is a rod-shaped bacterium, 0.7 µm in diameter, with a length usually in the range 2–6 µm. The swimming speed within a nominally uniform culture varies widely, but is typically less than 60 µm s⁻¹, with a mean of approximately 20 µm s⁻¹ (Kessler and Wojciechowski, 1997). The rate of oxygen consumption is estimated to be 10⁻⁶ molecules cell⁻¹ s⁻¹ (Berg, 1983). Oxygen uptake lowers the concentration in the main body of the culture, but diffusion from the air/water interface generates a concentration gradient. Individual bacteria sense the oxygen concentration gradient, swim up it and, therefore, accumulate near the top of the fluid. This accumulation results in an increase in the local average density of the suspension, since the cell density of B. subtilis exceeds the water density by approximately 10 %; Hart and Edwards (1987) give a value of 1.117±0.004 g ml⁻¹. In this situation, bioconvection is thought to begin as a gravitational Rayleigh–Taylor instability (Sharp, 1984), and its progression can be described as a sequence of states (Jánosi et al., 1998). (i) The well-mixed initial state. (ii) A horizontally uniform vertical oxygen concentration gradient develops just below the interface between the suspension and the air above it. Bacteria located within the gradient region start to swim upwards. (iii) Bacteria accumulate near the air/water interface, leaving a region depleted of bacteria just below. There is no visible instability (i.e. horizontal variations in cell concentration). (iv) In the fluid near the interface, the concentrated, relatively high mass density layer of bacteria-charged water ‘wrinkles’ because of the gravitational instability. Descending regions contain many cells; because of incompressibility, they ‘pull’ down the surrounding upper layers of fluid, thereby creating a positive feedback. The fluid adjacent to these downwelling regions also moves downwards because of the viscous stress. Incompressibility requires the generation of upwelling regions interspersed between the descending ones. (v) The plumes descend, and the extra cells drawn or swimming into them increase the opacity in bright-field visualization or increase the scattering (whiteness) in dark-field visualization. Cells tend to swim into the downwelling regions since they contain oxygen advected from the air/water interface region. (vi) The plumes hit the bottom and spread, like a typical gravity current; upwelling fluid compensates for descending plumes. (vii) Steady bioconvection is approached. The wavelength changes from that characterizing the ‘onset’, i.e. gravitational instability, to one that is appropriate for steady bioconvection (Bees and Hill, 1997). The geometry/symmetry also changes accordingly.

States i–iii do not produce changes observable from above; state v yields maximum contrast; state vi leading to state vii characterizes the onset of the steady nonlinear regime that develops from the spread of the plumes’ feet to generate rolling convection. Until the beginning of state iv, the leading anisotropy, i.e. the oxygen gradient, is vertical, as is the orienting gravitational field in the case of gravitaxis-mediated bioconvection.

One goal of the present study was to provide systematic data pertaining to bioconvection patterns due to consumption and upward chemotaxis for comparison with future theoretical predictions and numerical simulations. Our work is a natural extension of that of Bees and Hill (1997), who performed similar studies on the gravitactic bioconvection of the alga Chlamydomonas nivalis. Thus, another goal was to reveal systematic differences between these two varieties of bioconvection.

Despite the seemingly generic nature of the hydrodynamic instability, definite differences were found in the control parameter dependence and time evolution of gravitactic and aerotactic bioconvection patterns. In the case of gravitaxis, the characteristic scale of the convection cells is approximately equal to the depth of the suspension, as Fig. 7 demonstrates. In contrast, for chemotactic bioconvection, we found that it is the cell density that determines the dominant wavelength, while suspension depth plays a less important role. To quantify thisdifference, we fitted the two-parameter linear regression derived in the present study (equation 6) to the data of Bees and Hill (1997) data, yielding

What is the source of the differences in the coefficients A and B between the present results and those of Bees and Hill (1997) that distinguish algal and bacterial bioconvection? In the algal case, all the cells, throughout the entire culture, swim upwards because of gravitaxis. In the bacterial case, only those cells that sense the oxygen gradient (see above) swim upwards. When the cells that are near the air interface swim upwards, they leave behind a region of fluid depleted of cells (see above). The effective depth of the dense and depleted layers causing the onset of the instability is usually much shallower than the entire depth of the fluid, and the gradient in mass density is enhanced by the contrast between accretion and depletion. The phenomenon, unobservable in plan view, is easily seen from the side, as illustrated by Fig. 15.3 in Kessler and Wojciechowski (1997).

The reported depth-dependence of bacterial bioconvection patterns in tilted dishes [see Fig. 1 in Hillesdon et al. (1995); Fig. 1 in Metcalfe and Pedley (1998)] does not contradict our observations (Fig. 5B). In our analysis, the patterns at the onset of bioconvection were evaluated (see Jánosi et al., 1998), corresponding to region B in Fig. 3h. In tilted dishes, longer wavelengths in the deeper part of the tilted container were observed in the stationary regime (state vii; see above), corresponding to regions D–G in Fig. 3h.

The time evolution of the convection patterns is also markedly different for algae and bacteria. The patterns become more fine-scaled with time in the case of C. nivalis, equivalent to a decreasing wavelength (see Figs 2 and 8 in Bees and Hill, 1997). For B. subtilis, the opposite process occurs (Fig. 3). This results from the penetration of the pattern from the shallow accumulation/depletion region near the air interface to the entire depth of the culture. This penetration of convection cells entails a change in symmetry, compared with the onset of bioconvection: initially, the gradient of oxygen concentration is vertical, but it eventually becomes skewed away from the vertical by entrainment in the cellular convection. Conversely, one may infer that the trend to finer scales in the gravitactic case may result either from a net upward accumulation of cells, thereby producing a shallower convecting region, or from the generation of ‘bottom-standing plumes’ by gyrotaxis, swimming oriented by gravity together with vorticity (Kessler, 1985; Pedley and Kessler, 1992). Bottom-standing plume patterns are the eventual steady state in deep layers that can easily be observed from the side. It is not possible to obtain side views of shallow layers.

What about the shape of the wavelength spectra at and near the onset of bioconvection? Cursory inspection of the photographs (Fig. 1A, Fig. 3b and especially the magnified Fig. 6a) reveals a plethora of scales. The smallest is due to the video display and its interference with the main image. This has the right dimensions to account for the points between 1.0 and 2.0 on the abscissa of the plot in Fig. 4, which should be ignored. The narrowest ‘real’ scale is represented by the bright lines that are the curtains at the core of the plumes that descend and spread along the bottom. The principal additional scales are the separations between curtains, the diameters of the spreading plumes and the septa spacings between the plumes. Many additional experiments will have to be performed to unravel the relationship between the wavelength spectrum and these scales. This task will be accomplished using swimming mutants and other species of bacteria whose swimming styles produce variously shaped patterns (J. O. Kessler, unpublished observations) as well as suitably controlled initial conditions.

We note that the spectrum I(k) is closely related to the structure factor, which is basically the Fourier-transform of the pair-correlation function (equation 2). The structure factor is usually measured in the context of scattering experiments from physical objects. Disordered structures, random aggregates or fractal objects typically exhibit structure factors with power-law decays (Sinha, 1989; Hamburger-Lidar, 1996). Although our evaluation of images of photons scattered from a patterned array of organisms does not correspond to the usual scattering experiments, it may still be plausible to expect an intensity spectrum with a power-law decay.

The initial time-dependence of λ₀ (from a to b) shown in Fig. 3h demonstrates the natural development of sharp features out of an initially uniformly distributed, low-optical-contrast, unstable blanket of high mass density. The actual beginning of the instability will be situated somewhere along the extrapolation of the initial points to earlier times.

Future experimental studies will investigate the dependence of the time evolution and geometry of patterns on controls that relate biology, e.g. the age distributions of cells and their swimming velocity probability density function, with changes in the physics, e.g. temperature and viscosity, and the nutritional/biochemical background, e.g. growth medium constituents.

This work was supported by NATO Collaborative Research Grant No. CRG-960634, the Hungarian National Science Foundation (OTKA) under grants No. F026645 and T032423, and US National Science Foundation grant DMR-9812526. J.O.K. also gratefully acknowledges the support of Alice and Ralph Sheets, via a grant to the University of Arizona Foundation.