## Abstract

In a one-variable, finite size reaction–diffusion system, the existence of a minimal domain size required for the existence of a non-zero steady state is predicted, provided that the reaction–diffusion variable has a fixed value of zero at the boundaries of the domain (Dirichlet boundary conditions). This type of reaction diffusion model can be applied in population biology, in which the finite domain of the system represents a refuge where individuals can live normally immersed in a desert, or region where the conditions are so unfavourable that individuals cannot live in it. Building on a suggestion by Kenkre and Kuperman, and using non-chemotactic *E. coli* populations and a quasi-one-dimensional experimental design, we were able to find a minimal size (approximately 0.8 cm) for a refuge immersed in a region irradiated with intense UV light. The observed minimal size is in reasonable agreement with theory.

## 1. Introduction

*Escherichia coli* and other bacterial species have been observed to form a great variety of patterns under various environmental conditions and stresses (Ben-Jacob *et al*. 1995, 2000; Budrene & Berg 1995; Matsushita *et al*. 1998). In an effort to understand the collective dynamics of bacteria, various authors have proposed mathematical models in which reaction–diffusion is usually the model of choice, as reaction–diffusion models are capable of pattern formation through the Turing instability. In this paper, we present experimental confirmation of a theorem for a general reaction–diffusion equation in one-variable using populations of a non-chemotactic strain of *E. coli* bacteria. Building on a suggestion by Kenkre & Kuperman (2003), in which the authors proposed a way to test the validity of the Fisher equation (a special case of the reaction–diffusion model we considered) when applied to bacterial populations, we used a slightly modified version of an experimental setup initially used by Lin *et al*. (2004) and show that our results can verify the general reaction–diffusion theorem.

A simple one-variable reaction–diffusion system consists of a single equation,(1.1)where *c*(*x*, *t*) is a local concentration. The first term on the right corresponds to the random motion of particles or individuals, where *D* is the coefficient of linear diffusion. The second term corresponds to the reaction function *f*(*c*). When applied to population biology, *c*(*x*, *t*) corresponds to a single species population and *f*(*c*) corresponds to the growth (and death) of the population. In the case of a finite domain in which *c*(*x*, *t*) has a fixed value of 0 at the boundaries of the domain (Dirichlet boundary conditions with value 0), the finite domain is interpreted as a refuge, in which a population can grow normally, immersed in a desert where the population cannot live and an individual dies as soon as it leaves the refuge. Given these boundary conditions in a one-dimensional space, and if *f*(*c*) meets the conditions *f*(0)=0 and *f*′(0)>01 it has been shown that the refuge size should be greater than a critical size in order to support a non-null steady state (Skellam 1951; Murray 1993; Kenkre & Kuperman 2003) or, in biological terms, a steady population. The value of the critical domain size depends on the diffusion and growth function of the system, and is given by(1.2)where .

In any realistic model, the growth function should make the population saturate, as the population cannot grow forever in a closed system. One common saturating growth function is logistic growth,(1.3)where *a* corresponds to the linear birth rate and *K* corresponds to the maximum concentration possible in the system, called the carrying capacity. In this special case, the reaction–diffusion equation is known as the Fisher equation, which was introduced to model the spread of an advantageous gene in a population (Fisher 1937). Skellam (1951) and Kenkre & Kuperman (2003) showed how it is possible to find exact solutions of the steady state of the Fisher equation with Dirichlet boundary conditions in terms of Jacobi elliptic functions, and, from these solutions, the same result for the size of a critical domain (equation (1.2)) is obtained.

However, Dirichlet boundary conditions are an idealization of a real boundary, be it natural or fabricated. A more realistic approach to study critical domain sizes would be to use an equation in which individuals outside the refuge die at a certain rate so that the population outside the refuge need not be zero, but could account for individuals who are dying. Decay and diffusion could create a steady population outside the refuge. An example of such a system could be a modified Fisher equation in which the population grows logistically inside the refuge and dies at an exponential decay rate outside the refuge (we will show later that this is the case in our experiments). In this case, the modified Fisher equation would be(1.4)where *a* is the growth rate inside the refuge, *a*_{d} is the death rate outside the refuge, and *D* is the diffusion coefficient (assuming same *D* inside and outside the refuge).

In this case, no boundary conditions of the entire domain are specified, as we consider only the case where *a*_{d}>*a*. In this case, the bacteria die soon after leaving the refuge, and, therefore, they do not travel too far from the refuge boundary. If the system is much larger than the refuge itself, the boundary conditions do not affect the dynamics of the bacteria. It should also be noted that this simple model does not take into account the possibility of ‘damaged’ bacteria inside the refuge: some bacteria that go outside the refuge but later come back inside, before the UV irradiation kills them. This damaged bacteria may constitute a distinct population (with different traits) inside the refuge, for which the model does not account.

Ludwig *et al*. (1979) showed that for the model of equation (1.4) the size of a critical domain would become(1.5)This equation shows that if *a*_{d}→∞, then this equation becomes equation (1.2). If *a*_{d} has the same value of *a*, the size of the critical domain is reduced by half. Note that the same result applies to any function *f*(*c*) that meets the conditions *f*(0)=0 and *f*′(0)>0, since the result is obtained by linear stability analysis.

To investigate the existence of a critical domain size in a biological system we used an experimental setup initially used by Lin *et al*. (2004), in which the refuge was used to verify the existence of a phase transition between localization and extinction in *E. coli* populations when there was a convective flux that drove the bacteria outside the refuge.

In this setup bacterial populations live in a quasi-one-dimensional channel that is irradiated from the top with intense UV light, except for a small region in the centre, which is shielded from the UV light with a mask. We used *E. coli* bacteria for the experiments because of the fast growth of *E. coli* in microbiological media, which allowed us to conduct an experiment in a reasonable amount of time. Since the bacteria do not die instantly under the UV light, the conditions are more similar to equation (1.4) than to Dirichlet boundary conditions.

In this work, we compared the experimental results with equation (1.5). In §2, we present a detailed description of the experimental setup. In §3, we show how we obtained the parameters *a*, *a*_{d}, *D* and *K* experimentally in independent measurements. In §4, we present the results obtained, including the patterns obtained for various refuge sizes; in §5, we discuss the significance of the results and their implications.

## 2. Materials and methods

The experimental setup consisted of a glass channel, constructed with two glass pieces of 1.25×10×1/16 in. that sandwiched a Plexiglas spacer in a U shape, creating a rectangular container of dimensions 23×1.2×0.3 cm. The channel was held upright by a metallic base. Under the metallic base, a device that held a laser diode aligned with a photodetector could be moved along the direction of the channel to obtain concentration profiles in space. A schematic is presented in figure 1.

All this setup was placed inside an incubator to keep the temperature at 37 °C. On top of the incubator, and right above the channel, a mercury UV lamp (Ster-L-ray GT486L, *λ*=254 nm) was placed to illuminate the channel. The UV light was able to penetrate the incubator through a fused silica window. To achieve a collimated, uniform illumination of the channel, a black lighting control grid (10° angle spread) was placed between the fused silica window and the channel. The intensity produced by the UV lamp at 5 cm (the distance between the top of the channel and the lamp) from the lamp was 1750±80 μW cm^{−2}. However, the collimator reduced the intensity to 23% of the original, so the UV intensity at the top of the channel was 400 μW cm^{−2}. The channel was built with the length much greater than its width and height to create a quasi-one-dimensional space so that we could compare concentration profiles with one-dimensional models. The channel was filled with 5 ml (so that the actual height of liquid was 0.72 cm) of M9 medium supplemented with 4 g l^{−1} casamino acids, 1 mg l^{−1} Thiamine and 0.1% agar. M9 was chosen instead of other media because it is colourless, and UV light can penetrate it. With LB broth or tryptone broth, the UV light penetration is very poor. M9 supplemented with casamino acids and vitamins allowed fast growth, almost as fast as with nutrient broths. The agar was added to the solution to avoid external convection while still allowing bacteria to swim. To avoid evaporation of the medium, the incubator was held at high relative humidity (more than 75%), and a thin layer of mineral oil was applied above the medium before inoculation.

To obtain concentration profiles (concentration versus distance), the laser/photodetector scanned the channel length in 2500 steps, 0.095 mm apart, in a total time of 8 min and 30 s. In each step, the amount of light from the laser beam transmitted through the bacterial suspension was measured, the same way that a spectrophotometer works. The logarithm of the ratio of the initial and transmitted light is proportional to the concentration at that point, according to the Beer–Lambert law. The signal from the photodetector was fed into a computer to produce concentration versus distance data. To ascertain that the bacterial concentration was linear with the photodetector signal in the whole range of concentrations, we calibrated the photodetector signal using known dilutions and direct plate counts. The calibration is shown in figure 2. Our laser wavelength was 670 nm. We note that the concentration was linear with the signal over the whole range of concentrations in our case (as opposed to linearity in a spectrophotometer), because the width of the channel was 0.3 cm, as opposed to 1 cm width in standard cuvets.

The motion of the bacteria *E. coli* has been studied extensively (Berg & Brown 1972; Berg 2000). In a medium free of attractants, bacteria randomly alternate ‘runs’ (quasi-straightforward motion) with ‘tumbles’ (random change of direction without translation), making their overall motion a random walk. However, in the presence of attractants (as in most microbiological media), the motion of bacteria is biased toward attractant molecules. Therefore, the view of the motion of bacteria as a random walk is not clear anymore. The above considerations justify our choice of *E. coli* strains.

At first we wanted to use the *E. coli* strain RW120 (Ho *et al*. 1993), the same strain used by Lin *et al*. (2004), because the DNA repair is defective in this strain. Thus, the death rate under the UV lamp would be higher than normal. However, RW120 is a wild-type of *E. coli* for chemotaxis (this means that even though the genome has been modified for DNA-repair mechanisms, its mechanisms for chemotaxis remain the same as naturally occurring *E. coli*), and after growing for some time in the centre of the channel, *E. coli* RW120 forms a couple of chemotactic bands. These bands are swarms of bacteria that move outward very quickly, in response to their attraction to serine and aspartate. Adler (1966) was the first to observe and analyse such chemotactic bands; Keller & Segel (1971) were the first to propose a mathematical model for these bands. The model of Keller and Segel was later shown to be the continuous limit description of the motion of individuals in a reinforced random walk, similar to the trajectories observed for chemotactic *E. coli* under a microscope. (See Hillen & Othmer (2000) and references therein.)

In molecular biology terms, chemotaxis is understood to be a consequence of *E. coli* chemoreceptor proteins binding to attractant molecules (Parkinson 1993; Berg 2000). Mathematically, chemotactic bands correspond to travelling pulses, which cannot be solutions to equation (1.1). An equation like equation (1.1) allows the formation of travelling fronts, as the population spreads out, but not travelling pulses (Edelstein-Keshet *et al*. 1998). The observation of travelling pulses (figure 3) means that the bacteria do not move randomly and, therefore, the system could not be a simple one-variable reaction–diffusion system. We then decided to use the strain RP9535 (Sanatinia *et al*. 1995) which cannot respond to chemotactic signals, and whose motion is free from tumbles: every change in direction is due to collisions with other cells or molecules. This makes the motion of bacteria more random, making the reaction–diffusion model more valid, even though the death rate under the UV light is about eight times lower than for RW120, and the velocity of spread from the centre in our medium is about 1/10 of that of RW120. See figure 3.

## 3. Parameter estimation

To calculate the diffusion coefficient of RP9535, 10 ml of saturated solution were spun down in a microcentrifuge to reduce the total volume to 0.5 ml. A volume of 20 μl of this highly concentrated solution was slowly pipetted down at one end of the glass channel, which was filled with fresh medium. We observed the spread of the bacteria from the border for only 1.5 h, to minimize the effect of bacterial growth. By calculating the mean-square-distance of spread from the end point of the initial high concentration at the end of the channel versus time and comparing this mean-square-distance versus time with numerical simulations of the diffusion equation , we estimated *D*=(2.2±0.15)×10^{−5} cm^{2} s^{−1}. Other estimates of *E. coli* diffusion coefficients were done by Lin *et al*. (2004), although in that case the estimation was done for *E. coli* RW120, whose motion can be far from random, as explained above.

To calculate the value of the death rate under the UV light, as well as the carrying capacity of the medium, we used a glass channel with the same width and height of the original channel, but whose length was only 6 cm. We placed 2 ml of saturated solution into this small channel and then turned on the UV light. Every 30 min, a sample of 10 μl was taken after stirring the solution to assure uniform concentration. The concentration was calculated by the standard method of plate counting. From the graph of bacterial concentration versus time, the death rate was *a*_{d}=(9.46±0.45)×10^{−4} s^{−1}. The carrying capacity of the medium can be simply estimated as the maximum bacterial concentration observed in a test tube, which was *K*=(5.5±0.4)×10^{8} cells ml^{−1}. Note that the carrying capacity has no effect on the critical size calculation.

The growth rate of bacteria was measured from the total growth observed in the experimental setup, with the UV light switched off. The total growth was the sum of the concentrations at each spatial point. From the growth curve obtained, the growth rate was *a*=(2.23±0.2)×10^{−4} s^{−1}. The linear fits used to obtain the parameters are shown in figure 4.

## 4. Results

Each experiment was begun by dipping a sterile needle into a saturated bacterial solution and then submerging the tip of the needle into the centre of the channel, which was filled with microbiological medium (see §2 for details). The UV light illuminated the channel uniformly from the top, except for a small region in the middle where a mask blocked the UV light (the refuge). The size of masks that we used were 0.8, 1, 1.2, 1.5, 2 and 4 cm. We did three runs for each mask. The growth and dynamics of the bacteria were monitored for a period of 48 to more than 100 h. The bacterial population grew and spread from the inoculation point. At each side of the mask the bacteria moved outside the shaded region and kept spreading until the UV light killed the population outside the refuge. The growth curve of the total population (inside and outside the refuge) is similar to that of a test tube in which, after a time delay, the bacteria grow exponentially for a few hours, but then eventually saturates and remains constant for several more hours before the population starts to decay. We monitored each experiment during the time required for the concentration pattern to be approximately stable. This time was reached several hours after the total population saturated, since the pattern kept changing after the total population remained constant. To monitor the stability of the pattern, we used the *Χ*^{2} comparison test between one profile and the previous profile. When the *Χ*^{2}-test remained close to zero for several hours after the total population stopped growing, the experiment was stopped. Typical concentration profile dynamics and characteristics for the strain RP9535 are shown in figure 5.

Using the parameters obtained in §3 and equation (1.5), the value of the critical mask size was *L*_{c}=0.71±0.08 cm. Under the UV light intensity that we had, the bacteria went 1.5–2 cm away from the border of the mask before stopping. The growth under the mask invariably produced an alternating pattern, which is not reproduced exactly in every run. Each pattern had approximately 2–3 peaks cm^{−1}. The alternating patterns in the profiles were visible to the naked eye by looking very closely at the glass channel after an experiment was done, as the bacteria aggregated in vertical stripes under the mask. For the 2 and 4 cm masks, the time before growth showed in the signal varied between 3 and 5 h, a time comparable to that for bacterial growth in homogeneous space. For smaller mask sizes this time became larger, as a critical mask size was being approached. For the 1.5 cm mask the average time before growth was 12 h; for the 1.2 cm the average time was 25 h. In this last case, there was large variation between trials (10, 16 and 49 h). The 1 cm mask was very close to the critical size, and, therefore, the delay times were expected to increase significantly. One of the 1 cm mask's trials did not grow at all after 70 h; in another trial, the bacteria took 80 h to begin growing and grew to a maximum concentration whitin the range of maximum concentrations seen for other patterns (10^{9}–1.2×10^{9} cells ml^{−1}). The maximum concentration for each mask size was approximately constant. This finding disagrees with the Fisher model, which predicts that the maximum concentration decreases as the mask size decreases (Kekre & Kuperman 2003). The critical mask size was experimentally found to be 0.8 cm<*L*_{c}<1 cm, since no growth was observed in the 0.8 cm mask for a time of 110 h, during three trials. This value is 11% higher than 0.71 cm, the value obtained theoretically. The general characteristics of the patterns observed are shown in figure 6.

As a control experiment, we performed an experiment with *E. coli* strain RW120, a wild-type for chemotaxis, under the 0.8 cm mask. In this case, *E. coli* RW120 grew under the 0.8 cm mask after a time delay of 32 h. See figure 7. In this case, in addition to some other differences, the quasi-static pattern seems to be less stable than with the other strain. The distance travelled outside the mask was similar to RP9535, but the population that went outside the refuge was significantly smaller. Moreover, it was observed that the two peaks of maximum concentration inside the refuge are about 2 mm from the border of the mask as opposed to the RP9535 peaks, which invariably were right on the boundaries of the refuge. The fact that RW120 grows inside the 0.8 cm mask emphasizes the different behaviour of the two *E. coli* strains and confirms that equations (1.1) and (1.4) are not a reasonable model for chemotactically driven organisms.

## 5. Discussion

We have observed the existence of a minimal refuge size for the sustained growth of non-chemotactic *E. coli*. The minimum refuge size corresponds to the minimal domain in which the bacteria can have sustained growth and reach a non-zero, quasi-stable concentration profile. The value obtained experimentally was 0.8±0.1 cm, in contrast to the 0.71±0.08 cm obtained theoretically. *Escherichia coli* RW120, which has normal chemotaxis mechanisms, was observed forming chemotactic bands in homogeneous space (no refuge/desert situation) and growing in the 0.8 cm refuge, indicating that a one-variable reaction–diffusion equation is not a reasonable model for chemotactic bacteria. Consequently, the minimal refuge size theorem cannot be verified with chemotactic *E. coli*.

At first glance, it may seem difficult to imagine why the bacteria *E. coli* RP9535 would not grow in a 0.8 cm refuge where there is space to grow, as the RW120 strain shows.

To explain this, it is better to consider the microscopic situation rather than the reaction–diffusion model. Each RP9535 bacterium in Brownian motion has an average net displacement proportional to the square root of time (). This average displacement corresponds to the average distance travelled from the centre (in any direction) by a bacterium in a certain amount of time. As time passes, each bacterium has a larger average net displacement; Consequently the probability that it remains inside the refuge decreases. If there were no growth inside the refuge, all bacteria would eventually leave the refuge and die. Therefore, the situation can be viewed as a competition between growth and displacement. To have sustained growth inside the refuge, the growth must outbalance the displacement of bacteria toward the outside of the refuge. To compare growth and displacement in our experiments, we considered the average residence time in the refuge (, where Δ*x* is half the size of the refuge) and the average doubling time of a bacterium (, where *a* is the growth rate). With *a*=2.23×10^{−4} s^{−1}, the doubling time is *T*_{D}=51.8 min. As for the average residence time, we have *T*_{R}=94.7 min in the 1 cm refuge, and *T*_{R}=60.6 min in the 0.8 cm refuge. In the 0.8 cm refuge, the residence time is very close to the doubling time, meaning that most bacteria do not have time to reproduce before leaving the refuge. Furthermore, the doubling time as calculated corresponds to the lower limit of the doubling time, because *a* is calculated when the bacteria grow exponentially, which is only for a few hours, as shown in figure 6*a*. After that, the growth rate decreases, and the doubling time increases. The residence time can then become significantly less than the doubling time. Murray (1993) and Kenkre & Kuperman (2004) postulate that the critical refuge size in the reaction–diffusion model is a consequence of the competition between growth and loss at the boundaries; this explanation is consistent with the model of a finite number of random walkers that reproduce only inside the refuge.

As for the strain RW120 (which is a wild-type for chemotaxis) growing inside the 0.8 cm refuge, we note that the explanation above does not apply as they do not move randomly; a bacterium may not necessarily tend to leave the refuge. The fact that these bacteria grew in the 0.8 cm refuge is also consistent with the fact that both in our experiments (figure 3) and in Adler's original experiments in capillary tubes (Adler 1966), a significant number of bacteria remained around the inoculation point (in a non-motile state) after the chemotactic bands had formed and travelled outward.

The fact that the critical domain size theorem is verified with non-chemotactic *E. coli* does not imply that the one-variable reaction–diffusion model (equations (1.1) and (1.4)) is an accurate model for RP9535 growth at every stage. As shown in figure 3, the initial growth of RP9535 seems consistent with diffusion and growth only. When the total growth is reduced, the bacteria tend to form multi-peaked patterns in space. To model pattern formation more complex models would be needed. The one-variable reaction–diffusion model could be an approximation for the initial stages of growth. This is consistent with the fact that the pattern formation of bacterial colonies seems to be associated with lack of nutrients (Ben-Jacob *et al*. 1995, 2000; Matsushita *et al*. 1998). At the initial stage of growth there are plenty of nutrients available, and the bacteria reproduce fast. As the nutrients become scarce, pattern formation is observed as the bacteria compete for resources and employ survival mechanisms.

A one-variable reaction–diffusion system (equation (1.1)) describes what is essentially a system of non-interacting particles: each bacterium moves around randomly and reproduces according to *f*(*c*), independently of each other. However, chemotactic bacteria communicate with each other. This can be either indirectly, via the local concentration of chemoattractants initially present in the medium, or directly since under certain circumstances *E. coli* can excrete chemoattractants themselves (Budrene & Berg 1995; Mittal *et al*. 2003; Park *et al*. 2003). As a consequence, bacteria tend to form aggregates of high concentration. We found that a critical size is present for non-chemotactic *E. coli*, which implied that non-chemotactic bacteria behave like non-interacting particles (at least when there are plenty of nutrients available). On the other hand, the result cannot be verified with chemotactic bacteria, which suggest that they interact with each other. The idea that chemotaxis is the main way by which *E. coli* bacteria communicate with each other is supported by Park *et al*. (2003), and the main result from this paper reinforces that idea.

In summary, we have shown that the prediction of a critical domain size in one-variable reaction–diffusion systems can be verified in non-chemotactic *E. coli* populations. This implies that non-chemotactic *E. coli* behave as non-interacting particles, as opposed to *E. coli* with normal chemotaxis mechanisms.

## Acknowledgments

This work was carried out in the Laboratory of Chemical and Biological Patterns at Duke University. The author thanks Maggie Whitson and Grete Tornow for their advice in *E. coli* culturing, Daniele Armaleo for advice, John S. Parkinson for kindly providing the strain RP9535 and for his advice, Eric Monson for technical help and helpful discussions and Philippe M. Binder for his help in writing this document. This work was supported in part by Research Corporation (Award CC5885).

## Footnotes

↵ Author and address for correspondence: Biophysics Research Division, University of Michigan. 930 N University Avenue, Ann Arbor, MI 48109-1055, USA

↵These conditions are standard for population growth without an external source, and state that a population can only grow from a non-zero initial population.

- Received March 8, 2005.
- Accepted June 1, 2005.

- © 2005 The Royal Society