Theoretical Investigation of Bacterial Chemotaxis in Porous Media

We also thank Rajamani Rajagopalan, Jeannie Darby, Ed Schroeder, Brian Wood, and the ...... Jang-Kyoo Shin , Dong-Sun Kim , Hey-Jung Park , Geunbae Li...
0 downloads 0 Views 214KB Size
Download PDF
5636

Langmuir 2001, 17, 5636-5645

Theoretical Investigation of Bacterial Chemotaxis in Porous Media Kirk E. Nelson* and Timothy R. Ginn Department of Civil and Environmental Engineering, University of California, Davis, California 95616 Received March 26, 2001. In Final Form: June 18, 2001 Chemotaxis is defined as the self-generated displacement of motile microbes in the direction of an increasing concentration gradient of a chemoattractant or a decreasing concentration gradient of a chemorepellent. This work presents a Monte Carlo computer simulation of the chemotactic response in porous media. The model is based on the idealized Happel cell model underlying the colloid filtration theory (CFT), which describes the physical and chemical forces acting on an inert particle in flow through porous media. The premise of the model presented here is that bacterial chemotaxis in porous media can be successfully modeled by integrating cellular dynamics simulations with CFT particle trajectory analysis. At each time step a stochastic algorithm chooses a chemotactic velocity to be superimposed onto the deterministic velocity obtained by analytical solution for Stokes flow around the Happel sphere. The model calculates a value for the collection efficiency (the fundamental CFT parameter η) with the effects of chemotaxis incorporated. The collection efficiency is defined as the fraction of particles approaching the Happel sphere, termed the collector, that actually collides with the collector. CFT has been used extensively for the modeling of particle deposition in water treatment processes and, more recently, bacterial transport in groundwater. However, there has been little rigorous examination of the theory since it was developed in the mid-1970s. Thus, in addition to the treatment of chemotaxis, this paper provides an analysis of the theoretical underpinnings of CFT. Preliminary results suggest that chemotaxis may significantly affect η under certain conditions.

I. Introduction A. Filtration Theory. Predicting the fate and transport of microorganisms in porous media is an important problem because of the health hazards posed by pathogenic bacteria and viruses and also because of the potentially beneficial use of microbes in the bioremediation (in particular, biostimulation and bioaugmentation) of contaminated aquifers. The modeling of particle transport (and removal) in filter media is a complex problem even neglecting biological processes such as chemotaxis. Over the past few decades, the body of work known as colloid filtration theory (CFT) has been developed. For an accessible and detailed review of CFT, readers are referred to Rajagopalan and Tien.1 Despite its many simplifying assumptions (see Table 1), CFT has proved to be a useful modeling tool in practice. To apply CFT to microbial transport in groundwater, one may incorporate the collection efficiency η into a transport model as a reaction rate kinetic

{

∂M ∂2M ∂M 31- ) DL 2 - vx + (Rη) ∂t ∂x 2 d ∂x

[

]}

(1)

where M is microbe concentration, t is time, DL is the dispersion coefficient, x is distance, vx is the average linear groundwater velocity,  is porosity, d is the average diameter of porous media grain, R is the collision efficiency (defined as the fraction of particles colliding with the collector that actually attaches to its surface), and η is the collection efficiency (defined as the fraction of particles approaching the collector that actually collides with the collector). * To whom correspondence should be addressed. Tel: (530) 7529677. Fax: (530) 752-7872. E-mail: [email protected]. (1) Rajagopalan, R.; Tien, C. Prog. Filtr. Sep. 1979, 1, 179-269.

Initially, CFT was used to describe the rate of colloidal particle deposition during vertical-packed bed filtration in water treatment. The use of CFT to model bacterial transport in groundwater aquifers was first suggested2 in 1988 and implemented3,4 in the early 1990s. The initial approach has been to calculate η a priori and then to use this value in the calibration of the collision efficiency, R, via column experiments. The consequent values of R and η are then inserted into an expression for an adsorption constant. Of relevance here is that η is calculated by formulas resulting from a deterministic analysis. For example, a force and torque balance was performed on the suspended particle yielding a differential equation, the solution of which gives η.5 The first step of the current work was to calculate the clean-bed collection efficiency (so named because the effects of particle buildup are neglected) of CFT in a way compatible with a stochastic analysis of chemotaxis. Rajagopalan and Tien’s calculation of η is based on the trajectory concept (see Table 1 for a list of some of the assumptions of this model).5 With Brownian motion neglected, the remaining forces are deterministic: interception, sedimentation, and surface forces. Thus, for any starting position on the Happel shell, an exact trajectory can be computed when stochastic parts of the motion such as Brownian displacement are neglected. It is not necessary to compute all possible trajectories (which is impossible anyway since there exist an infinite number of possible trajectories) to determine η. Because η is defined geometrically as the ratio of the circular cross section of flow encompassing all particles (e.g., bacteria) that collide (2) Mathess, G.; Pekdeger, A. J. Contam. Hydrol. 1988, 2, 171-188. (3) Harvey, R. W.; Garabedian, S. P. Environ. Sci. Technol. 1991, 25, 178-185. (4) Martin, R. E.; Bouwer, E. J.; Hanna, L. M. Environ. Sci. Technol. 1992, 26 (5), 1053-1058. (5) Rajagopalan, R.; Tien, C. AIChE J. 1976, 22, 523-533.

10.1021/la010456+ CCC: $20.00 © 2001 American Chemical Society Published on Web 08/04/2001

Bacterial Chemotaxis in Porous Media

Langmuir, Vol. 17, No. 18, 2001 5637

Table 1. Notes on the Use of Rajagopalan and Tien5 for a Sub-Pore-Scale Analysis of Colloid Transport aspect of R & T (1976) model

significance

flow field around grains is assumed to be accurately represented by Happel model

since actual flow is random and non-uniform at and below the pore scale, the model can be considered an approximation of average flow at best if initial particle deposition has a significant effect on subsequent deposition, the model is invalid if ap < 1 µm, deterministic analysis must be combined with stochastic analysis; equations 18 and 27 in Rajagopalan and Tien give an inconsistent description of interception, sedimentation, and diffusion; a corrected expression exists in the literaturea model is not conceptually consistent with nonuniform groundwater flow unless neutral buoyancy of bacteria is assumed natural media grains are neither smooth nor spherical, most microbes would be better represented as ellipsoids

the collector surface is assumed to remain clean Brownian diffusion is assumed negligible if particle radius is >1 µm flow is assumed to be vertical and downward particle is collected when separation distance equals zero, both particle and collector are assumed to be perfectly smooth and spherical collector is at least two orders of magnitude larger than particle, allowing collector surface to be assumed flat double layer force is neglected in the approximate closed form solution for η the expression for K2 in the stream function is missing a minus sign the y variable is incorporated in the angular component of drag and torque correction factors aRajagopalan,

London attraction term and drag correction factors are invalid if condition for the assumption is not met the surface force is assumed to be solely determined by the London-van der Waals force use of stream function as written will result in large errors the order of y must be decreased/increased by one when writing uθ with/without drag and torque correction factors

R.; Tien, C.; Pfeffer, R.; Tardos, G. AIChE J. 1982, 28, 871-872.

Figure 1. For interception alone, the limiting trajectory ends at (as + ap, π/2). When other forces (e.g., London-van der Waals force, Brownian motion) are considered, particles may collect at angles >π/2. 2-D spherical coordinates (i.e., polar) represent a 3-D flow field under the assumption of axisymmetric flow.

with the collector to the circular cross section of total flow passing through the Happel shell, one only needs to find the so-called limiting trajectory, which separates the intercepting trajectories from those that pass on by. Then the collection efficiency is simply calculated from the initial angular coordinate θs (see Figure 1) of the limiting trajectory

η ) sin2(θs)

(2)

Because chemotaxis is treated as a stochastic component of the motion, the addition of a chemotactic velocity negates the utility of the limiting trajectory concept. Instead, it is necessary to use a stochastic approach, such as to compute a sufficiently large number of trajectories iteratively in Monte Carlo fashion. Such an approach is presented here and shown in the absence of chemotaxis to be consistent

with the results of Rajagopalan and Tien5 for interception and the London-van der Waals attraction force under the conditions tested. B. Chemotaxis. Bacterial chemotaxis is the process by which motile bacteria move toward more favorable conditions in response to chemical gradients. Motile cells, propelled by their flagella, swim in the direction of their long axis at a rate of about 35 diameters/s, often changing course but rarely stopping.6 Chemotaxis was poorly understood until the groundbreaking work of Berg and Brown7 with Escherichia coli. With the use of a microscope that automatically follows individual cells,8 it was discovered that a cell’s motion consists of an alternating sequence of “runs”, in which the cell swims at a nearly uniform velocity, and “twiddles” (later called “tumbles” by other researchers), in which the cell makes an abrupt change of direction according to a probability distribution skewed toward small angles.7 In the absence of a chemical gradient, the cell prefers no direction and the result is random motility. It may be possible that random motility in porous media can be modeled with an effective random motility coefficient (analogous to the effective diffusion coefficient used for the transition regime between bulk and Knudsen diffusion).9 However, it should be noted that bacterial motility differs from molecular diffusion in the nonuniform distribution of turn angles and the finite time required to make a turn. In the presence of a chemical gradient, Berg and Brown observed that the distribution of turn angles remained constant regardless of the chemical gradient being experienced by the chemotactic strains of E. coli. However, the probability of terminating a run decreased when a cell was moving in a favorable direction and increased when a cell was moving in an unfavorable direction. Thus, it is through the adjustment of run times, and not turn angles, that the bacterium biases its movement. This probabilistic structure results in a dispersive transport that is not the same as Brownian motion with drift, which supposes an additive drift (6) Berg, H. C. Phys. Today 2000, 24-29. (7) Berg, H. C.; Brown, D. A. Nature 1972, 239, 500-504. (8) Berg, H. C. Rev. Sci. Instrum. 1971, 42, 868. (9) Duffy, K. J.; Cummings, P. T.; Ford, R. M. Biophys. J. 1995, 68, 800-806.

5638

Langmuir, Vol. 17, No. 18, 2001

displacement that effectively changes both the angle and duration of a displacement step. Several mathematical models have been formulated to describe bacterial chemotaxis in a bulk aqueous medium.10-14 Of these, the most general and rigorous is Alt’s model of three-dimensional (3-D) cell balance equations.13 Unfortunately, the complexity of this model precludes an analytic solution. A thorough treatise on the relationships between these various models has been presented in the literature.15 This treatise also derives a reduced form of the Alt model under the assumption of symmetry in two of the three coordinate directions. One way to solve this reduced Alt model is to simulate the motion of a large number of cells.16 This approach has been termed cellular dynamics simulation (CDS) by way of analogy to molecular dynamics simulation. Duffy et al. modified the CDS algorithm for a bulk aqueous medium to simulate the random motility of Pseudomonas putida (a bacterial species known to degrade certain groundwater contaminants) in the presence of an idealized porous medium.9 The simulations began with 1000 bacterial cells placed uniformly throughout the pore space. If the displacement over a time step placed a cell within the solid fraction of the porous medium, its original position was maintained. The simulation results showed that tortuosity (i.e., the ratio of random motility in the bulk aqueous medium to the effective value in the porous medium) increased as obstacle diameter decreased. This trend translates into cell motility becoming increasingly short-circuited as the porous matrix is constructed of smaller and smaller grains. Duffy et al. used CDS to model chemotaxis in porous media by adding a radial chemoattractant gradient to the random motility simulations and replacing the constant tumbling probability with a concentration dependent probability adapted from the one used for a bulk aqueous medium.17 The expression for the tumbling probability includes a parameter called the chemotactic sensitivity coefficient (χo), which is a measure of the bacterial flux caused by the chemotactic response. The higher χο is, the longer a cell will run when it is heading in the direction of an increasing attractant gradient. The simulations suggested that chemotaxis can enhance bacterial migration through porous media as long as the population’s χo remains below a certain threshold value. A large value of χo, which would enhance the chemotactic response in a bulk aqueous medium, may in fact short-circuit chemotaxis in a porous medium. The reason for this paradox is that longer runs may result in increased collisions with solid surfaces. While impeding transport, this condition could help bacteria increase their residence time in the vicinity of a food source. So it is possible that there exists both a threshold χo, not to be exceeded if enhanced transport is sought, and a minimum χo that must be met to increase residence time in the vicinity of a food source. The model presented below may be useful in evaluating the effects of χo, because the Happel cell lends itself to the analysis of local concentration gradients within pores. (10) Keller, E. F.; Segel, L. A. J.Theor. Biol. 1971, 30, 225-234. (11) Stroock, D. W. Z. Wahr. verw. 1974, 28, 305-315. (12) Segel, L. SIAM J. Appl. Math. 1977, 32, 653-665. (13) Alt, W. J. Math. Biol. 1980, 9, 147-177. (14) Rivero, M. A.; Tranquillo, R. T.; Buettner, H. M.; Lauffenburger, D. A. Chem. Eng. Sci. 1989, 44 (12), 2881-2897. (15) Ford, R. M.; Cummings, P. T. SIAM J. Appl. Math. 1992, 52 (5), 1426-1441. (16) Frymier, P. D.; Ford R. M.; Cummings, P. T. Chem. Eng. Sci. 1993, 48, 687-699. (17) Duffy, K. J.; Ford, R. M.; Cummings, P. T. Biophys. J. 1997, 73, 2930-2936.

Nelson and Ginn

II. Model Development The simulations reported thus far in the literature17 indicate some important relationships that may exist between the chemotactic behavior of bacteria and the geometry of a porous soil matrix. However, the effects of convective flow and the physical-chemical interactions that occur between porous media surfaces and bacterial surfaces have not been considered. CFT has modeled both of these factors to some degree of success.3,4 Therefore, the premise of the model presented here is that bacterial chemotaxis in porous media can be successfully modeled by integrating CDS with CFT particle trajectory analysis. In sections A through E, we treat the deterministic components of CFT (i.e., everything except Brownian motion): interception, London-van der Waals attraction, hydrodynamic retardation effect, double layer force, and sedimentation. The Monte Carlo approach to calculating η is developed in section F. Finally, the algorithm for chemotaxis is presented in section G. A. Interception. Collection, due to interception, results from the finite size of particles. Particles are assumed to follow streamlines exactly. If a particle is riding on a streamline that comes within one particle radius of the collector, then it will be collected. The streamlines are found by transforming the stream function of the Happel sphere-in-cell porous media model18 (which describes the spherical collector moving upward through a stationary fluid) to the case in which the collector is stationary and the liquid is moving down.5 Assuming incompressible, axisymmetric creeping flow, the resulting stream function is

ψ(r,θ) )

[

()

]

K1 U + K2r* + K3r*2 + K4r*4 as2 sin2(θ) 2 r* (3)

where ψ is the stream function, r is the radial coordinate, θ is the angular coordinate, U is the approach velocity of the liquid, as is the collector radius, r* is the dimensionless radial coordinate (r* ) r/as), and K1, K2, K3, and K4 are functions of porosity. The Darcy flux should be used for U. It should be noted that the expression for K2 (see Nomenclature) is missing a minus sign in the original paper deriving this part of CFT.5 The formula for K2 has appeared correctly19 after the orginal mistake, but to the writers’ knowledge the existence of the error in the original paper has not been pointed out previously in the literature. (Also worth noting is that the literature contains a clarification of aspects of CFT that are relevant for the macroscopic use of the theory.20 See Table 1 for discussion and clarification of aspects that are relevant to a subpore-scale analysis.) Taking the appropriate partial derivatives of eq 3 and dividing through by r sin θ (since the area of flow is proportional to r sin θ), the velocity of liquid around the collector is given by

[

v(r,θ) ) -

] [

]

1 1 ∂ψ ∂ψ er + e ∂θ r sin θ ∂r θ r sin θ 2

(4)

where er is the unit vector in r direction and eθ is the unit (18) Happel, J. AIChE J. 1958, 4 (2), 197-201. (19) Tien, C. Granular Filtration of Aerosols and Hydrosols; Butterworth: Boston, 1989. (20) Logan, B. E.; Jewett, D. G.; Arnold, R. G.; Bouwer, E. J.; O’Melia, C. R. J. Environ. Eng. (N.Y.) 1995, 121, 869-873.

Bacterial Chemotaxis in Porous Media

Langmuir, Vol. 17, No. 18, 2001 5639

For the particle tracking algorithm, it is necessary to have these expressions in their dimensional form as follows

ur ) -A+(1 + δ+)2U uθ ) [B+(1 + δ+) + D+(1 + δ+)2]

(8) U r

(9)

A detailed account of how the A+, B+, and D+ coefficients are obtained is not given in the literature. Thus, the following procedure was developed. First, an expression for A+ can be found algebraically from eqs 4 and 5

A+ ) Figure 2. Liquid velocity of the Happel model approaches zero as θ f 0 and as θ f π because pore throat diameter f ∞ at these angular coordinates, termed “stationary points” of the flow field. Velocity approaches its maximum as θ f π/2 because pore throat diameter reaches its minimum at π/2.

vector in θ direction. Figure 2 gives a surface plot of the velocity field around the collector. Later in the CFT model development, drag and torque correction factors are introduced to account for deviations from Stokes law near the collector surface.5 These correction factors require eq 4 to be re-expressed in the approximate form

v ) -Ay2er + [By + Dy2] eθ

(5)

where y is the distance from the particle centroid to the collector surface. It should be noted that the angular component drag and torque correction factors (not presented in this paper) incorporate the y variable. This explains why y is reduced by 1 order in the CFT angular particle velocity expression of Rajagopalan and Tien.5 A, B, and D are coefficients obtained by approximating eq 4 in the form of eq 5 over segments of the liquid shell. Thus, A, B, and D are functions of (r,θ) but are expressed as piecewise constants over segments of the liquid shell. The resulting particle velocity expressions for interception alone (omitting the drag correction factors) are

1 dr ) -A+(1 + δ+)2 u r* ) U dt

( )( ) r dθ u * ) ( )( ) ) B (1 + δ ) + D (1 + δ ) U dt +

θ

+

+

+ 2

(6) (7)

where ur * is the dimensionless radial component of velocity, uθ* is the dimensionless angular component of velocity, A+ is dimensionless A, B+ is dimensionless B, D+ is dimensionless D, and δ+ is the dimensionless separation distance between particle and collector. In the particle tracking algorithm described below, the expression for the normalized separation distance (δ+) is not allowed to be less than 0.001, corresponding to the bacterium being one one-thousandth of its own radius away from the collector. It is assumed that bacteria spend little time in this interval and that other complicating forces (e.g., DLVO) at this separation do not come into play. However, this assumption may be tenuous in light of recent results showing measurable electrosteric interactions over distances on the order of bacterial radii.21 (21) Camesano, T. A.; Logan, B. E. Environ. Sci. Technol. 2000, 34, 3354-3362.

{(

)( )( )}

2 1 ∂ψ ap y2r2 sin θ ∂θ U

(10)

A+ is undefined at θ ) 0 and θ ) π. To accommodate this behavior, θ is not allowed to be smaller than 0.001 rad. It does not appear tractable to determine B+ and D+ analytically, because there exists only one equation (By + Dy2 ) (1/r sin θ)(∂ψ/∂r)) for these two unknowns. Thus, the following approximation scheme was developed. For any r coordinate inside the liquid envelope, a segment is defined as having radial length 2∆y such that

b - as Ny

∆y )

(11)

where Ny equals twice the total number of segments. In the initial runs of the particle tracking program, Ny was set at a value equal to 104. For a particle at (r,θ) inside the liquid envelope, the angular component of eq 4 is first evaluated at (r - ∆y,θ) with the conditions that ∆y is set to zero if y is less than or equal to ∆y and θ is set to 0.001 rad if θ ) 0 or π. Then it is evaluated at (r + ∆y, θ). The first computed velocity and the y value associated with it are called v1 and y1. The second ones are v2 and y2. Linear interpolation yields

D+ )

B+ )

( )

( ) 2

ap U

( )(

v2 v1 y2 y1 y2 - y1

)

ap ν2 - Dy22 U y2

(12)

(13)

Equations 3, 4, and 8-13 define the particle velocities for interception. We now briefly describe the testing of this portion of the model, the results of which appear in Figure 3. The collection efficiency due to interception alone can be obtained directly from the stream function.5 This is the fractional approach area of streamline start locations on which colloids of radius ap will by virtue of their size (assumed spherical) strike the collector. The exact expression in terms of the stream function can be closely approximated by

ηI )

() ( )

as 2 3 NR2 As 2 b

(14)

where ηI is the collection efficiency due to interception alone, NR is the relative size group (NR ) ap/as), and As is a parameter accounting for the effect of neighboring grains (a function of ). Calculations of η using the particle tracking model were compared to values obtained from eq 14. The parameters

5640

Langmuir, Vol. 17, No. 18, 2001

Nelson and Ginn

that arises from the force and torque balance is

NLO )

H 9πµap2U

(15)

where H is the Hamaker constant and µ is the fluid viscosity. This expression is accurate only for extremely small separation distances. When the particles are separated by a distance equal to or greater than λe, the wavelength of electron oscillation (∼0.1 µm), there is a delay in the interactions between molecules causing a decrease in the attractive force. Thus, NLO must be scaled by R(r), the retardation factor, which accounts for the decrease in the attractive force as separation distance increases. An expression for R(r) was developed for use in CFT. Readers are referred to the literature for the full expression.22 Note that this R is not to be confused with the collision efficiency. The London attraction term is added only to the particle’s radial velocity expression, because in the idealized Happel sphere model the London-van der Waals force acts along the normal direction and, thus, has no effect on angular velocity. The radial velocity expression for combined interception and London attraction is

[

ur-LO ) -A+(1 + δ+)2 -

Figure 3. Interception collection efficiency as a function of porosity, particle radius, and collector radius.

that determine collection due to interception are particle radius, collector radius, and porosity. Values of each parameter were varied with the other two held constant. The model was found to be consistent with the results of CFT. The results are given in Figure 3 and Table 2. B. London Attraction. The London-van der Waals force is caused by momentary pairwise oscillations in electron density when two particles are near each other. In each particle, the electron clouds surrounding atomic nuclei are repelled causing electron density to be momentarily greater in one region of the particle than in another. This condition creates instantaneous dipoles that draw the particles together. The expression used in CFT to account for the Londonvan der Waals force5 was derived under the assumption that the spherical collector can be considered as a flat surface.19 This assumption is considered valid if the filter grain is at least 2 orders of magnitude larger than the approaching particle (i.e., 0.01 g NR, an assumption that may be violated for some natural porous media with distributed grain size). The dimensionless London group

R(r)NLO +2

]

δ (2 + δ+)2

U

(16)

The set of expressions for interception and London attraction was tested for  ) 0.39, as ) 100 µm, U ) 1 m/day, and NLO ) 1.45 × 10-5 with ap set equal to 0.1, 0.59, and 1 µm (corresponding to NR ) 10-3, 5.9 × 10-3, and 10-2). A comparison of the results obtained here with those given by Rajagopalan and Tien5 in their Figure 2 for NG ) 0 and drag corrections excluded is presented in Table 3 and Figure 4. The slight increase in collection (compared to the simulations of interception alone) is consistent with the intuitive expectation of the effects of an attractive force. C. Hydrodynamic Retardation Effect. The drag and torque correction factors that appear in the trajectory equation5 of CFT account for the hydrodynamic retardation effect. This effect is caused by the presence of a solid boundary near the suspended particle. The result of the effect is an increased viscous resistance to particle transport near the collector that is not accounted for in the streamline model. The literature contains both tabulated values23 and analytical expressions24-28 for these factors. The analytical expressions involve infinite series terms that are not easily truncated. Therefore, the tabulated values were fitted to develop expressions for this model. The results obtained are compared to the values produced by eq 25 of Rajagopalan and Tien5 in Table 3 and Figure 4. The drag and torque correction factors were not used in the chemotaxis simulations. D. Double Layer Force. The double layer force is a surface force, the sign of which depends on the signs of (22) Payatakes, A. C. A new model for granular porous media: application to filtration through packed beds. Ph.D. Dissertation, Syracuse University, Syracuse, NY, 1973. (23) Rajagopalan, R. Stochastic modeling and experimental analysis of particle transport in water filtration. Ph.D. Dissertation, Syracuse University, Syracuse, NY, 1974. (24) Brenner, H. Chem. Eng. Sci. 1961, 16, 242-251. (25) O’Neill, M. E. Matematika 1964, 11, 67. (26) Goldman, A. J.; Cox, R. G.; Brenner, H. Chem. Eng. Sci. 1967, 22, 637-651. (27) Goldman, A. J., Cox, R. G.; Brenner, H. Chem. Eng. Sci. 1967, 22, 653-660. (28) Goren, S. L. O’Neill, M. E. Chem. Eng. Sci. 1971, 26, 325-338.

Bacterial Chemotaxis in Porous Media

Langmuir, Vol. 17, No. 18, 2001 5641

Table 2. Interception Collection Efficiency as a Function of Particle Radius, Collector Radius, and Porosity  ) 0.39, as ) 100 µm

 ) 0.39, ap ) 0.59 µm

ap ) 0.59 µm, as ) 100 µm

ap (µm)

ηmodel

ηR&T76

as (µm)

ηmodel

ηR&T76



ηmodel

ηR&T76

0.1 0.2 0.3 0.4 0.5 0.59 0.6 0.7 0.8 0.9 1.0

4.32 × 10-5 1.69 × 10-4 4.00 × 10-4 6.76 × 10-4 1.09 × 10-3 1.52 × 10-3 1.52 × 10-3 2.12 × 10-3 2.70 × 10-3 3.48 × 10-3 4.22 × 10-3

4.36 × 10-5 1.74 × 10-4 3.93 × 10-4 6.98 × 10-4 1.09 × 10-3 1.52 × 10-3 1.57 × 10-3 2.14 × 10-3 2.79 × 10-3 3.53 × 10-3 4.36 × 10-3

5 15 25 45 65 85 100 150 250 350 500

0.493 0.063 0.023 7.38 × 10-3 3.60 × 10-3 2.11 × 10-3 1.52 × 10-3 6.76 × 10-4 2.56 × 10-4 1.21 × 10-4 6.03 × 10-5

0.607 0.067 0.024 7.50 × 10-3 3.59 × 10-3 2.10 × 10-3 1.52 × 10-3 6.75 × 10-4 2.43 × 10-4 1.24 × 10-4 6.07 × 10-5

0.20 0.23 0.26 0.31 0.37 0.41 0.43 0.49 0.51 0.53 0.60

8.08 × 10-3 5.92 × 10-3 4.35 × 10-3 2.81 × 10-3 1.76 × 10-3 1.30 × 10-3 1.09 × 10-3 7.29 × 10-4 6.76 × 10-4 5.76 × 10-4 3.61 × 10-4

8.46 × 10-3 6.05 × 10-3 4.48 × 10-3 2.85 × 10-3 1.76 × 10-3 1.31 × 10-3 1.14 × 10-3 7.55 × 10-4 6.61 × 10-4 5.79 × 10-4 3.67 × 10-4

Table 3. Collection Due to Interception, London Attraction, and Hydrodynamic Retardation Effect without drag & torque corr’ns

with drag & torque corr’ns

NR

ηR&T76

ηmodel

ηR&T76

ηmodel

0.001 0.0059 0.010

5.62 × 10-5 1.38 × 10-3 4.07 × 10-3

6.40 × 10-5 1.76 × 10-3 4.75 × 10-3

1.70 × 10-5 4.80 × 10-4 1.30 × 10-3

4.90 × 10-5 7.80 × 10-4 1.20 × 10-3

the particle and collector surface potentials. All solid materials placed in an aqueous medium acquire a surface charge, which is balanced by countercharged ions present in the solution. This double layer of charge gives rise to an electrical potential between the outer portion of the double layer and the bulk solution, the zeta potential, which can be used to approximate the potential difference between the material surface and the bulk solution. Recent results29 suggest that this approach may be erroneous. The double layer force may be either repulsive or attractive depending on the signs of the particle and collector surface potentials. The double layer force and the London-van der Waals force comprise the total potential described by DLVO theory (see Ninham, for a review of the capabilities and limitations of DLVO theory30). The assumption that surface interactions are dictated by DLVO theory may not be accurate for bacterial colloids.21 Rajagopalan and Tien5 neglected the double layer force in their approximate closed form solution for η, and it has been neglected in this model as well. E. Sedimentation. Although it appears straightforward to incorporate the effects of gravity into the particle tracking algorithm, it was not done for two major reasons. First, as previously pointed out,4 the orientation of the Happel sphere-in-cell porous media model makes it conceptually consistent only with vertical flow. However, we are often interested in horizontal (or generally nonuniform) groundwater flow for bioremediation schemes and natural colloid transport modeling. If the effects of gravity are neglected, there is no difference in collision rate between horizontal and vertical flow. Second, while the specific gravity of native and laboratory-cultured bacteria is a matter of some debate, it is generally accepted that indigenous bacteria are less dense than their cultured counterparts. Without a reliable way to measure indigenous bacterial densities, it seems reasonable to assume neutral buoyancy. Under this assumption, collection due to sedimentation would be zero. For these reasons, sedimentation was neglected. However, a sensitivity analysis of density could indicate if there is reason to incorporate sedimentation into the model. (29) Elimelech, M.; Masahiko, N.; Ko, C. H.; Ryan, J. N. Environ. Sci. Technol. 2000, 34, 2143-2148. (30) Ninham, B. W. Adv. Colloid Interface Sci. 1999, 83, 1-17.

Figure 4. Collection due to interception and London attraction, with and without the hydrodynamic retardation effect.

F. Calculating η without the Limiting Trajectory Concept. The random effects of chemotaxis invalidate the limiting trajectory concept. It is possible that a bacterium with a greater starting angle will collect and one with a lesser starting angle will not. Therefore, the Monte Carlo approach described below was developed. The range of possible starting angles is from 0 to π/2 rad (see Figure 1). A very small positive angle called θo is chosen for the first trajectory to be computed. Trajectories are computed from n different starting angles between θo and π/2. The Happel sphere-in-cell is actually a 3-D model represented with two-dimensional (2-D) polar coordinates under the assumption of axisymmetric flow. Thus, each starting angle, θi, represents not a single point in 2-D space but rather a circular cross section of fluid. (The starting angles are projected onto the approaching tube of fluid above the collector so that they represent circular cross sections comprised wholly of fluid. However, the circular cross sections corresponding to the flow field around the collector are comprised of both fluid and the solid material of the collector.) Assuming a uniform particle concentration, the sine of each θi gives the number density of particles at that angle. To make density weights that are all integers, the first angle is chosen by

θo ) sin-1

(n1)

(17)

All of the starting angles are then given by

θi ) sin-1 (i sin θo)

(18)

i ) 1, 2, ..., n And the density weights are

Ni )

sin θi sin θo

(19)

5642

Langmuir, Vol. 17, No. 18, 2001

Nelson and Ginn

Figure 5. Flowchart of simulation algorithm.

This scheme results in N1 ) 1, N2 ) 2, and so on. Neglecting chemotaxis, one trajectory can be run for each θi, and if it collects, all particles starting from that angle are considered collected. Computing η in this manner for  ) 0.37, NR ) 0.012, NLO ) 1.45 × 10-5, and n ) 100, a value of 3.77 × 10-3 was obtained. When the limiting trajectory was found to three decimal place accuracy and the resulting θs was used in eq 2, an η of 7.039 × 10-3 was obtained. This discrepancy indicates the level of truncation error. G. Chemotaxis. Under the assumption that the chemotactic velocity and the velocity resulting from the forces described by CFT are additive, a separate expression for the bacterium’s velocity due to chemotaxis was developed. At each time step, (∆t), the program first checks if a tumble occurs. If a tumble occurs, a new swimming direction is chosen and the bacterium is transported under CFT alone for the duration of the tumble. Then, the chemotactic and CFT expressions are added to obtain an overall velocity from which the displacement over ∆t is computed. At the end of each time step, it is checked to see if the bacterium has left the Happel cell (i.e., an uncollected trajectory) or reached the collector surface (i.e., a collected trajectory). A separate trajectory is run for each of the Ni bacteria starting from a given θi. After every trajectory is run, the process is repeated for a desired number of realizations. The computation of chemotactic velocities is based on the Duffy et al. simulations17 of P. putida. Figure 5 gives a conceptual flowchart of the entire algorithm for calculating the chemotactic η. The decision to turn is based on the concentration dependent turn frequency (eq 8 in Duffy et al.17)

[

β(r,θ) ) βo exp

]

Kd -χo s‚∇C(r,θ) v (K + C(r,θ))2 d

(20)

where βo is the turn frequency in the absence of a chemical gradient, χo is chemotactic sensitivity, v is the constant cell swimming speed, Kd is the dissociation constant for bacterial receptor-attractant binding, C is the chemoattractant concentration, and s is the unit direction vector in polar coordinates.

Figure 6. Chemotactic direction vectors.

In the simulations of Duffy et al.,17 bacterial motion is modeled as being entirely self-generated, equating s to the chemotactic direction. In this model, motion is caused by the combination of swimming and the physical and chemical forces described by CFT. Therefore, s in this case is the overall direction in which the bacterium is moving as a result of swimming, advection, etc. Note that the likelihood of a turn decreases as the s component of ∇C increases. The radial and angular components, respectively, of the chemotactic direction vector are

sr ) -cos(θ + γ)

(21)

sθ ) sin(θ + γ)

(22)

where γ is the last angle in which the bacterium has turned (see Figure 6). The overall expressions for radial and angular velocity are

vr ) ur-LO + (srv) vθ ) uθ +

(sθv) r

The overall unit direction vector, then, is

(23) (24)

[ ]

Bacterial Chemotaxis in Porous Media

Langmuir, Vol. 17, No. 18, 2001 5643

vr

s)

xvr2 + (vθr)2

(25)

vθr

xvr2 + (vθr)2

The decision of whether to turn is made at each time step as follows. A random number (called F) is chosen on the interval [0,1]. A turn will be made if β∆t > F. Otherwise, the bacterium’s existing swimming direction is maintained. If a turn is made, its magnitude γ is chosen from the bacterium’s turn angle distribution. From the data points in a plot of a reported distribution for P. putida,31 we derived the following seventh degree polynomial by the method of least squares

Figure 7. Collection with and without chemotaxis as a function of flow velocity. Chemotaxis appears to change the negative dependence of η on groundwater velocity to a positive dependence under the simulated conditions. Table 4. Simulation Results

P(γ) ) -0.002 + 1.842γ - 2.2610γ2 + 1.180γ3 0.478γ4 + 0.200γ5 - 0.048γ6 + 0.0038γ7 (26) On a typical length scale in groundwater systems (e.g., U ) 1 m/day,  ) 0.37, as ) 100 µm), P. putida will usually make between zero and two turns in a trajectory through the Happel cell. Thus, the initial value given for γ will have a significant impact on the rate of deposition. First, simulations were run setting the initial value of γ equal to zero. This condition biases movement toward the collector, because the CFT trajectories curve around the collector. As expected, simulations run with this initial condition yielded values of η significantly greater than those obtained without chemotaxis. A second condition tested was to choose the initial value of γ such that the initial swimming direction of each bacterium is the same as the initial direction of its velocity computed from CFT. This condition also biases movement toward the collector but not as much as the first condition. The formula for this condition is

(

γi ) tan-1

uθ(b,θi)

)

-ur-LO(b,θi)

- θi

(27)

The model chemoattractant simulated was 3-chlorobenzoate (3CB). 3CB, also modeled by Duffy et al.,17 is a structural analogue to benzoate and uses the same receptor system to elicit a chemotactic response, but it is not metabolized by P. putida. Thus, secondary gradients created by consumption may be neglected. The concentration gradient was modeled as a linear function of r and independent of θ. Concentration increased in a linear manner from 0.65 mmol at the collector surface to 1 mmol at the edge of the liquid envelope. The following parameters were taken from Duffy et al.:17 βo ) 0.5 s-1, χo ) 2 × 10-4 cm2/s, Kd ) 0.5 × 10-3 mol/L, and v ) 44 µm/s. The duration of a tumble was 0.025 s. The number of starting angles was set equal to 25. From eq 19, it follows that a total of 325 bacteria were simulated.

Darcy flux, U

nonchemotactic η 10-2

0.56 m/day

1.80 ×

120 m/day

9.00 × 10-3

realizations

chemotactic η

1 10 25 1 10 25

6.00 × 10-3 9.50 × 10-3 1.20 × 10-2 1.50× 10-2 1.60× 10-2 1.60× 10-2

port32). For both cases, the CFT parameters were set as follows:  ) 0.37 and NR ) 0.012 (based on ap ) 0.59 µm and as ) 50 µm). The London group depends on U; the low U gives NLO ) 0.175, and the high U gives NLO ) 8.183 × 10-4. At the low U, runs were made for 1, 10, and 25 realizations. A realization is one possible outcome of a stochastic process (random function). In this case, the random function is the collection efficiency η. For one realization, 2 out of 325 trajectories were collected, yielding η ) 0.006. For 10 realizations, an average of 3.0875 out of 325 trajectories were collected, yielding η ) 0.0095. For 25 realizations, an average of 3.9 trajectories were collected, yielding η ) 0.012. These results are low relative to the value η ) 0.018 computed without chemotaxis. In fact, each set of realizations produced an η lower than the nonchemotactic case, suggesting that under these particular conditions chemotaxis may enhance bacterial transport. At the high U, runs of 1, 10, and 25 realizations yielded η values of 0.015, 0.016, and 0.016, respectively. These results are all greater than the value η ) 0.009 computed without chemotaxis. Conventional CFT predicts the decrease of η as groundwater velocity increases (as was computed by this model with chemotaxis neglected). However, the simulations of chemotaxis reported here show the opposite trend. These results are consistent with the hypothesis of Camesano and Logan32 that motility and chemotaxis may enhance bacterial transport more effectively at lower porewater velocities. The simulation results are summarized in Table 4 and Figure 7. IV. Discussion

III. Results To gain insight into the effect of groundwater velocity on the chemotactic response, simulations were run at two different velocities: U ) 0.56 m/day (in the upper range of groundwater velocities in coarse grained material) and U ) 120 m/day (an extremely high velocity but chosen for consistency with a previous study on the relationship between flow velocity, cell motility, and bacterial trans(31) Duffy, K. J.; Ford, R. M. J. Bacteriol. 1997, 179, 1428-1430.

We have presented a simulation model for calculating the CFT parameter η (collection efficiency) with the added effects of chemotaxis. The model was shown to be consistent with the results of Rajagopalan and Tien5 for the collection mechanisms of interception and the Londonvan der Waals force. With chemotaxis incorporated under a particular idealized chemoattractant gradient, the (32) Camesano, T. A.; Logan, B. E. Environ. Sci. Technol. 1998, 32, 1699-1708.

5644

Langmuir, Vol. 17, No. 18, 2001

simulations produced a lower value of η (relative to the nonchemotactic value) at a low porewater velocity and a higher value of η at a higher porewater velocity. These results suggest that chemotaxis may enhance bacterial transport at low groundwater velocities and can increase bacterial residence time in the vicinity of a contaminated zone and that the dependence of η on porewater velocity changes sign when chemotaxis is considered. The model presented in this paper may be a useful tool in understanding how chemotaxis can be a significant and possibly exploitable trait in subsurface bacteria, including those used for bioremediation. The simulations done thus far have been exploratory. The first step toward producing more meaningful results will be to conduct sensitivity analyses on the number of starting angles simulated and the number of realizations run. A second refinement worth investigating is the implementation of a reflecting boundary condition (i.e., if a bacterium leaves the liquid shell at a small angle, it may be more realistic to consider it to be entering the flow field of a neighboring collector). Also worth examining is how the hydrodynamic retardation effect, sedimentation, and Brownian diffusion affect results. Although the modeling of sedimentation will not be conceptually consistent for horizontal groundwater flow, it could be important for analyzing the transport characteristics of cultured bacteria with higher densities. If Brownian diffusion is deemed significant, it would be a second stochastic process in the model. Additionally, it might be interesting to perform simulations within a hyperboloidal constricted tube model33 (instead of the Happel model) to investigate the effects of pore blocking on chemotaxis. Another intriguing topic is the chemotactic sensitivity parameter, χo. A sensitivity analysis of χo could indicate the potential for enhanced transport to a contaminated zone versus the potential for increased residence time in a contaminated zone. Another possibility is that motility by itself (with or without a chemical gradient) plays the dominant role, and so a comparison of random motility simulations and chemotaxis simulations should be made. Since this model is based on the behavior of individual bacteria cells within an idealized pore space, it may be appropriate to coordinate the modeling effort with microscale experiments of bacteria in constructed porous media. Eventually, it is possible that the model could be be used as a screening tool for bioremediation schemes and as a generator of chemotactic η values for use in a macroscopic transport model. Nomenclature R ) collision efficiency ) fraction of particles colliding with collector that actually attaches to it R(r) ) retardation factor accounting for decrease in London-van der Waals force with distance ap ) particle radius as ) collector radius A, B, and D ) coefficients obtained by approximating eq 4 in the form of eq 5 over segments of the liquid shell of Happel model A+ ) dimensionless A ) A(ap2/U) As ) neighboring grains parameter ) [2(1 - p5)/w] b ) radius of Happel cell βo ) turn frequency in the absence of a chemical gradient β(r,θ) ) turn frequency in the presence of a chemical gradient B+ ) dimensionless B ) B(ap/U) C ) chemoattractant concentration (33) Venkatesan, M.; Rajagopalan, R. AIChE J. 1980, 26, 694-698.

Nelson and Ginn δ+ ) normalized surface-to-surface separation distance ) (y - ap)/ap d ) average diameter of porous media grain D+ ) dimensionless D ) D(ap2/U) DL ) dispersion coefficient  ) porosity er ) unit vector in r direction eθ ) unit vector in θ direction γ ) bacterial turn angle γi ) initial turn angle of bacteria at starting angle θi η ) collection efficiency ) fraction of particles approaching collector that actually collides with it ηI ) collection efficiency due to interception alone H ) Hamaker constant ≈ 10-13 ergs i ) index of starting angles K1 ) 1/w ) stream function coefficient K2 ) -(3 + 2p5)/w ) stream function coefficient K3 ) (2 + 3p5)/w ) stream function coefficient K4 ) -p5/w ) stream function coefficient Kd ) dissociation constant for bacterial receptor-attractant binding λe ) wavelength of electron oscillation ≈ 0.1 µm µ ) fluid viscosity ≈ 8.94 × 10-4 N s/m2 M ) microbe concentration n ) number of starting angles from which trajectories are computed Ni ) number of bacteria at starting angle θi NR ) relative size group ) ap/as NLO ) London group ) (H/9πµap2U) Ny/2 ) number of segments over which B+ and D+ are constant p ) variable in stream function coefficients ) 3(1 - )1/3 ) as/b P(γ) ) bacterial turn angle distribution ψ(r,θ) ) stream function θ ) angular coordinate θo ) smallest starting angle from which trajectories are computed θs ) initial angular coordinate of limiting trajectory r ) radial coordinate s ) unit direction vector of overall bacterial motion sr ) radial component of chemotactic direction vector sθ ) angular component of chemotactic direction vector t ) time ∆t ) time step U ) approach velocity of liquid (e.g., Darcy flux or specific discharge for groundwater) ur ) dimensional radial component of particle velocity for interception alone ur-LO ) dimensional radial component of particle velocity for interception and London attraction uθ ) dimensional angular component of particle velocity for interception alone u/r ) dimensionless radial component of particle velocity for interception alone u/θ ) dimensionless angular component of particle velocity for interception alone v(r,θ) ) liquid velocity field of Happel sphere-in-cell porous media model v ) bacterial swimming speed v1 ) v(r - ∆y,θ)‚eθ ) value of liquid angular velocity used to find B+ and D+ at (r,θ) v2 ) v(r + ∆y,θ)‚eθ ) value of liquid angular velocity used to find B+ and D+ at (r,θ) vr ) overall bacterial radial velocity vθ ) overall bacterial angular velocity vx ) average linear groundwater velocity w ) variable in stream function coefficients ) 2 - 3p + 3p5 - 2p6 χo ) chemotactic sensitivity coefficient x ) distance y ) distance from particle centroid to collector surface ) r - as

Bacterial Chemotaxis in Porous Media y1 ) r - ∆y - as ) value of y associated with v1 y2 ) r + ∆y - as ) value of y associated with v2 2∆y ) length of radial liquid shell segment over which B+ and D+ are constant

Acknowledgment. This research was supported by the U.S. Department of Energy, Natural and Accelerated Bioremediation Research Program, partly under the project entitled “The Influence of Heterogeneity and Growth on Microbial Transport in Saturated Porous

Langmuir, Vol. 17, No. 18, 2001 5645

Media” and partly under the project entitled “Data Analysis/Integration for Comparative Oxic/Anoxic Field Injections of Microbes”. The authors thank Yannis Dafalias for the use of his Apple Power Macintosh G3 and Jeannie Darby for the loan of Payatakes’ dissertation. We also thank Rajamani Rajagopalan, Jeannie Darby, Ed Schroeder, Brian Wood, and the anonymous reviewer for their valuable comments. LA010456+