Originally published as Biophys J. BioFAST on August 25, 2006.
doi:10.1529/biophysj.106.091314
Biophysical Journal 91:3640-3652 (2006)
© 2006 The Biophysical Society
Mesoscopic Modeling of Bacterial Flagellar Microhydrodynamics
Yeshitila Gebremichael,
Gary S. Ayton and
Gregory A. Voth
Center for Biophysical Modeling and Simulation, Department of Chemistry, University of Utah, Salt Lake City, Utah
Correspondence: Address reprint requests to Gregory A. Voth, Center for Biophysical Modeling and Simulation, Dept. of Chemistry, University of Utah, Salt Lake City, UT 84112. E-mail: voth{at}chem.utah.edu.
 |
ABSTRACT
|
|---|
A particle-based hybrid method of elastic network model and smooth-particle hydrodynamics has been employed to describe the propulsion of bacterial flagella in a viscous hydrodynamic environment. The method explicitly models the two aspects of bacterial propulsion that involve flagellar flexibility and long-range hydrodynamic interaction of low-Reynolds-number flow. The model further incorporates the molecular organization of the flagellar filament at a coarse-grained level in terms of the 11 protofilaments. Each of these protofilaments is represented by a collection of material points that represent the flagellin proteins. A computational model of a single flexible helical segment representing the filament of a bacterial flagellum is presented. The propulsive dynamics and the flow fields generated by the motion of the model filament are examined. The nature of flagellar deformation and the influence of hydrodynamics in determining the shape of deformations are examined based on the helical filament.
 |
INTRODUCTION
|
|---|
Motile bacteria such as Escherichia coli and Salmonella typhimurium are known to be sensitive to light, temperature, and chemical variations in their local aqueous environment (1
). The response to chemical stimuli, referred to as chemotaxis, is manifested by the bacteria swimming toward chemoattractants or away from chemorepellents. This is done by sensing the chemical gradient at a faster rate than the time it takes for direction changes due to Brownian motion (1
). An external organelle of locomotion known as the flagellum serves as a rotating propeller to aid in the "swimming" process. The flagellum consists of three parts: a reversible rotary motor, a hook, and a filament. The motor is a complex structure composed of proteins (1
,2
). It is located at the flagellum's anchor point that is embedded in the cell wall, beginning with the cytoplasm and ending at the outer membrane. It is driven by a transmembrane electric potential due to a proton gradient (for some bacteria that live at high pH, sodium ions are used instead) across the cell membrane (1
,2
). The motor can run in both clockwise and counterclockwise directions, generating a torque that rotates a helical filament attached to it via a short proximal hook. The hook functions as a flexible coupling or universal joint for the filament, transmitting the torque generated by the motor to the filament (see Namba and Vonderviszt (3
) for reviews).
The flagellar filament serves as a propeller by converting the rotary motion of the motor into a linear thrust (4
,5
). The filament is a self-assembling helical polymer constructed from a single protein flagellin. It is able to adopt different polymorphic conformations, including left-handed and right-handed supercoils (6
,7
). It has 11 monomers for every two turns of the one-start helix. Alternatively, the filament is conceived as 11 strands of protofilaments arranged around a hollow central channel (3
,8
,9
). These helical filaments are up to 15 µm long, whereas they are only 120250 Å in diameter (3
). Typical flagellar filaments of contour length from 10 to 15 µm are estimated to contain from 21,000 to 32,000 flagellin monomers, respectively (10
).
A bacterial cell may contain a single flagellum or multiple flagella attached to it at different sites of the cell body. For example, E. Coli typically contains up to six flagella per cell (2
). During chemotaxis, the swimming pattern of the bacteria alternates between "run" and "tumble" modes (11
). During the run mode the motor rotates counterclockwise (as it is viewed from outside the cell), and left-handed helical filaments bundle together (12
), providing a forward thrust to the cell. Periodically, the bundle falls apart when one or more of the flagella motors switch their direction of rotation, causing the cell to move erratically or tumble (13
). As a result, the propulsion of the cell is carried out through a serious of runs interrupted by tumbles (11
). When the tumbling event occurs, a complex process is initiated, where a chirality reversal that eventually transforms a left-handed helix into a right-handed helix is manifested (14
).
The propulsion of a bacterial flagellum involves the interplay of elasticity and hydrodynamics. The structure of the flagellum is dominated by strong elastic behavior, i.e., in terms of bending and other deformation modes. On the other hand, the surrounding solvent is governed by viscous interactions and the Newtonian constitutive relation. As such, modeling the composite flagellum/solvent system is a challenging task, and that is the focus of this article. In the past, a large number of analytical and numerical calculations have been carried out to describe the nature of bacterial propulsion and the flow surrounding it (15
29
). For example, many theoretical studies of the fluid mechanics of helical swimmers are aimed at deriving the behavior of helical propulsion from first principles (15
17
). These studies predict the motion of helical structures and a large body attached to them. Lighthill (18
) provided a mathematical description of a helical swimmer using slender-body theory (19
). In this theory, the flagellum is represented as a line distribution of singularities known as Stokeslets and dipoles. It takes advantage of the long slender shape of the filament and obtains an approximate solution for the flow around such bodies. Higdon (20
) combined analytical and numerical approaches to describe the motion of a helical flagellum with a spherical body. The flagellum was modeled using slender-body theory. Computationally, Ramia et al. (21
) employed the boundary-element method to study the motion of a helical flagellum with a spherical body. Kim et al. (22
) used slender-body theory to study the flow surrounding a helical filament.
In all of these works, the helical filaments are considered as rigid structures that rotate as a whole, or deform in a presumed wave form at all times. As a result, the elastic nature of the bacterial flagellum is ignored. Alternatively, in other studies (23
26
), the elastic nature of the filament is captured through the application of Kirchhoff's rod theory, whereas hydrodynamics is only approximated through the use of resistive-force theory (27
). In the resistive-force theory, the force on a short segment of the filament is calculated by considering the local drag force normal and tangential to the line element of the helical filament. In particular, the force per unit length is approximated by f = 
u
+
||u||, where u
and u|| are the normal and tangential components of the local rod velocity relative to the fluid, and 
and
|| are the transverse and the tangential resistance coefficients, respectively. In this regard, the resistive-force theory is a simple approach to use. However, it makes no prediction for the flow induced by the filament.
Furthermore, in most of these treatments the structure of the propeller appears only as helical lines or tubes (20
,22
,24
). A model that fully addresses the bacterial flagellum needs to account for the structural organization of the filament as well. This description is lacking in many of the models. An exception to this particular case is the coarse-grained continuum rod model of a bacterial flagellum that incorporates the structural organization of the filament in terms of the protofilaments (29
). Another recent development in this respect is the work of Cortez et al. (30
), based on the immersed boundary method. This model treats the flagellum as a helical hollow tube consisting of material points connected to each other by elastic springs. The material points serve as force-generating point sources that are advanced on the basis of a no-slip boundary condition. The method incorporates both the elastic and the hydrodynamic aspects of bacterial propulsion. The choice for the material points in this model is not typically representative of the 11 strands of the protofilaments. However, it captures the actual swimming of a bacterial flagellum through the coupling of flexibility and hydrodynamics.
In this work, the molecular organization of a bacterial flagellum is constructed in terms of the protofilaments that constitute the propeller. In particular, the filament is represented as a hollow helical tube of 11 strands of protofilaments. Each of these protofilaments is composed of material points that are regarded as particles representing flagellin proteins at a coarse-grained level. Then a model of a bacterial flagellum in a fluid is developed through the use of a hybrid method of the elastic network model (ENM) and the smooth-particle hydrodynamics (SPH). Many authors have utilized the ENM method to describe the global motion of large macromolecules such as proteins (31
). In these models the nodes of the network represent the amino acid residues and the linkers are springs that represent the interresidue potential stabilizing the folded conformations. At a macroscopic level, on the other hand, such a network has also been used in the immersed boundary method (30
,32
). Similarly, in the present model, each material point on the model filament is treated as a node and any neighboring points on the surface of the helical tube are connected by a network of springs with a stiffness constant S. The stiffness constant of the model is matched to the material property of the real flagellum by connecting the local property of the network with the overall elastic property of the flagellum. The overall elastic property can be found from experimental measurements of either the flexural rigidity A or the twist modulus C of the flagellum.
The flow around the model filament and the coupling between the surrounding fluid and the filament are studied through the use of SPH. This method is a fully Lagrangian method in which numerical solution is achieved without the use of a grid. The gridless nature of SPH provides an advantage over the traditional fluid dynamical methods that use finite-element or finite-difference methods in dealing with complex geometries or highly deformable boundaries. The SPH method was originally developed for astrophysical problems (33
,34
) and later applied to incompressible fluids through the use of equation of states that yield low Mach number M (35
).
The goal of this work is thus twofold: 1), to capture the physics of bacterial propulsion by developing a simplified model that accounts for both the elastic and the hydrodynamic aspects of the flagellum, as well as the structural organization of the filament; and 2), to investigate the flow field and the deformation behavior arising from the interaction of elasticity and hydrodynamics. The flow field generated by rotating helices has recently been studied using an experimental technique of particle image velocimetry (PIV) (36
). Theoretically, there are few studies, and a rigorous description of the flow field is lacking. The explicit hydrodynamic methodology employed in the model described here sheds light on the nature of flow surrounding the filament. The flow field is investigated under conditions where the model filament is propagating as in a real flagellum.
The article is organized as follows: the simulation is briefly described in the Methods section, where a description of the model filament, as well as of the hybrid ENM and SPH methods, is provided. Further details about the model can be found in Supplemental Material. In the Results section, the propulsive dynamics of bacterial flagellum as captured by the hybrid model is first presented. Then, the effect of flagellar rotation on the surrounding fluid is studied by examining the flow field generated by a swimming helical filament. The reverse effect, i.e., the effect of hydrodynamics on the deformation behavior of the flagellum, is also examined. In this way, the self-consistent dynamics of bacterial flagellum in a hydrodynamic environment is described. A conclusion is given in the last section.
 |
METHODS
|
|---|
Model
The bacterial flagellum is a hollow tube of helical filament that consists of 11 protofilaments. Each of these protofilaments is an assembly of flagellin proteins. The filament is modeled at a coarse-grained level by setting up a tubular structure of helical geometry that consists of material points representing the flagellin proteins. The length scale
of a particle corresponding to each material point is determined by noting that there are 11 flagellin molecules around the tube of the filament, each of which belongs to one of the 11 strands of protofilament. Each particle is represented by a sphere of diameter
=
D/11, where D is the diameter of the filament. For real flagella, D is typically in the range 1225 nm. Using D = 24 nm,
is estimated to be
= 6.85 nm. The mass mo of the particles is estimated from the molecular weight of flagellin. This value ranges between 45 and 53 kDa (37
). By choosing a molecular mass of 50 kDa, mo is estimated as mo = 8.31 x 1020 g. The timescale
is selected based on the smallest time needed for stable simulation of the composite filament/solvent system. Accordingly, for the range of speed of sound and viscosity used in our simulation, a stable simulation under hydrodynamic environment is found for a time step
t = 0.001
, where
= 109 s. In the article, all values are quoted in reduced units, i.e., length in units of
, mass in units of mo, and time in units of
.
The model described here consists of 40 particles on each protofilament. The total number of particles that make up the model filament is thus 440. This number is chosen to allow a computationally manageable number of fluid particles that cover the filament. With this number of particles, the aspect ratio L/a, where a is the radius of the filament, is smaller than that of the full-scale bacterial flagellum. This model may thus be conceived as a model of short flagellum. The filament prepared in this way is surrounded by 64,000 SPH particles that are placed on a regular cube of volume. For simplicity, the mass and size of the SPH particles have been chosen to be the same as the flagellum particles. The density of the filament is estimated as
1.14, whereas the fluid density is
1.006.
The geometry of the model filament is constructed by aligning a helical tube along the z axis. The position r(i) of any material point i is determined by
 | (1) |
where R and P are the radius of curvature and the pitch of the flagellar helix, respectively. a is the radius of the filament and 0
z(i)
L, where L = 40.0
is chosen as the axis length of the filament.
(i) is an equally spaced angle used to assign positions for the 11 material points around the surface of the tube. R and P are set as R = 1.0
and P = L/4, yielding a ratio R/P = 0.1. For a real flagellum, which can be found typically in any of the normal, semicoiled, or curly-1 polymorphic forms (1
,38
), R/P
0.087 for the normal E. Coli, where R
0.19 µm and P = 2.28 µm are estimated for this form (38
). The ratio R/P determines the pitch angle
by the relation
= tan1(2
R/P).
28.7 for the normal state, whereas
43.3 and 55 for the curly-1 and semicoiled polymorphic forms, respectively (see Table 1 of Kim et al. (39
) for the values of R and P). For the model described here,
= 32.1. The model is thus close to the normal polymorphic state.
It should be emphasized here that the length scale of the present model is much smaller than that of the real flagellum. This choice is made to allow for computational efficiency. To accommodate for a helical shape within the length scale of the model, the individual dimensions corresponding to the radius of curvature or the pitch of the helix are assigned values that differ from those of the real flagellum. However, as described above, the combination of these parameters yields a helical structure that is similar to the normal polymorphic state of flagellum, with the pitch angle serving as a reference quantity. More important, the elastic material behavior of the model flagellum is arranged such that the overall elastic property of the real flagellum measured in terms of flexural rigidity is attained. Therefore, the flagellum in our model may be perceived as a helical filament with material points organized to satisfy the overall elastic property of the real flagellum.
Simulation
Elastic network model
The filament is modeled based on the elastic network model. Each material point is treated as a node, and neighboring particles within a cutoff radius rc are connected by a network of springs. A single stiffness parameter is used for all particles. The use of a single stiffness parameter is justified, since the filament is comprised of identical protein molecules. The total potential energy E of all the particles in the filament is
 | (2) |
where Eij is the potential energy between any pair i and j. It is expressed in terms of a simple pairwise Hooke's law potential with a single value of S for all links,
 | (3) |
Here
ij is a step function that vanishes for rijo > rc. For our system, rc is chosen to be 2.0
. The stiffness constant, S, is parameterized based on the material property of the real flagellum, from which S = 10.17 x 103 dyn/cm is found (cf. Supplementary Material A). The force on i due to j is thus calculated as
 | (4) |
Smooth-particle hydrodynamics
The hydrodynamics component of the simulation is modeled using the SPH method. The SPH equation of motion for particle i corresponding to the Navier-Stokes momentum equation is expressed as (35
)
 | (5) |
where pi and pj are the pressures at particles positions ri and rj, respectively. Similarly,
i and
j are the mass densities of particles at positions ri and rj, respectively. W is the smooth-function kernel. For our model, W is expressed using the Lucy weight function (see Eq. B4, Supplementary Material B) with a cutoff distance h = 3
. The term
ij is an artificial viscosity first introduced to permit the modeling of shocks (40
). In general, this term produces a shear and bulk viscosity. It is expressed as
 | (6) |
Here,
is the nondimensional viscosity parameter, vij = vi vj, and rij = ri - rj. For low Mach number flow, ß is often set to zero (35
). The mass density
i of a particle at position ri can be evaluated using the relation given by Eq. B2 (Supplementary Material B). Here, an expression derived from the continuity equation (35
) is utilized, i.e.,
 | (7) |
This equation has the advantage of reducing the computational cost by allowing
to evolve concurrently with the particle velocities. Furthermore, the presence of spurious pressure gradients induced at a free surface can be avoided by using this expression (35
).
The equation of state used in the present model is based on the relation that was originally proposed by Batchelor (41
) and later applied to the traditional SPH model (35
),
 | (8) |
where
o is the initial density of the SPH particles. B and
are constants, where
= 7 is used as in Monaghan (35
). This choice causes the pressure to respond strongly to variations in density. B is chosen in such a way that the speed of sound is sufficiently large to ensure M
0.1 (cf. Supplementary Material B).
The swimming of a bacterial flagellum is within the regime of low Reynolds number, Re. The typical value of Re for bacterium in water is Re
104 or less, but there are also estimated values as high as Re
0.050.25 for the complex flagellar filament of Rhizobium lupini (see discussions in Trachtenberg et al. (42
)). In the simulation described here, the Reynolds number is set such that Re lies within these ranges, where Re is estimated from the relation Re
vl/
ch. Here, l is the characteristic length over which the flow varies, and v is the characteristic velocity of the flow (cf. Supplemental Material B). For the flow induced by the flagellar rotation, v is estimated based on the relation v =
l, where l = a, the radius of the filament. The interaction between the SPH fluid and the model filament is set up through a boundary condition that imposes a no-slip condition for the layer of fluid near the filament. This condition is approximated by applying a viscous interaction between the fluid and the filament (43
). Additionally, a pressure force arising from the density gradient near the boundary is exerted on the filament. More details about the development of the SPH method specific to this system can be found in Supplemental Material.
 |
RESULTS AND DISCUSSION
|
|---|
Propulsive dynamics
Bacterial dynamics is in the realm of low-Reynolds-number hydrodynamics. The remarkable aspect of motion at very low Reynolds number is that a reciprocal, or time-reversion invariant, motion generates no propulsion. According to the "scallop theorem", as Purcell nominally termed it, an organism that has only one degree of freedom, as in the motion of a scallop, or a rigid one-armed swimmer, cannot swim at low Re (44
). This is a consequence of the time reversibility of Stokes flow, in which the pattern of displacement is the same whether the organism moves fast or slowly during the power or recovery stroke. Motile microorganisms overcome this difficulty by breaking the time-reversal symmetry either through flexibility or cyclic motion. For example, eukaryotic flagella break this symmetry through oscillatory bending (45
,46
). Bacteria cells break the time-reversal symmetry through a cyclic motion of their flagella. This motion is generated by a continuous twist exerted at the base of the filament by the rotary motor. The twist force at the base is then transmitted to the rest of the filament through the elastic interaction. As the flagellum rotates, each piece of the filament experiences a viscous drag that sums up to yield a net thrust force that propels the cell. In this section, it will be shown how the hybrid model captures these essential physics of flagellar propulsion.
The propulsion of bacterial flagellum is studied through the model filament by applying a rotational torque at the base of the filament that spins the filament at a particular rotational speed
. This torque imitates the effect of motor rotation on the real flagellum. A typical rotational speed during bacterial chemotaxis is
= 100 Hz (38
,47
), but there is evidence that as high as 1700 Hz is possible (48
). Since our model incorporates the molecular motif of the flagellum at the level of flagellin, to observe any considerable propulsion, such rotational speeds demand significantly large CPU time. To speed up the simulation, higher rotational speeds are used. Nevertheless, the Reynolds number is adjusted to approximate the range found in real systems by selecting a proper speed of sound that compensates for the large
. The approach in setting Re as a target is conceptually similar to experiments done to study flagellar bundling using a macroscale model of bacterial flagella. In these experiments, the motors are typically rotated at a low speed of 0.1 Hz (39
) and 0.25 Hz (36
). However, the other parameters in the experiments are adjusted to yield a low Reynolds number with Re
103 (39
).
Different rotational frequencies f = 0.0034
1, 0.017
1, 0.034
1, 0.068
1, and 0.17
1 are used to investigate the propulsion behavior. The lowest of these values is three orders of magnitude higher than the largest motor rotation experimentally predicted. However, for each frequency, the speed of sound, c, that yields a low Mach number of M = 0.1 results in Re
0.006 for
= 10. These sets of parameters are used in the subsequent analysis. To specifically observe the propulsion, f = 0.17
1 is used. The other frequencies also give the same behavior, but takes longer time to cover the same displacement.
When the filament is rotated counterclockwise at a particular frequency
, the interaction between the filament and the fluid produces a thrust force that propels the flagellum along its axis, in a direction parallel to the direction of
. The filament moves in the opposite direction when the rotation is reversed. Fig. 1 shows a sequence of snapshots depicting the propulsion of model flagellum for the motor rotational frequency f = 0.17
1. The first frame shows the helical structure near the beginning of the simulation. To surround the model filament by fluid at all times, a periodic boundary condition has been applied along the long helical axis. The part of the filament that appears to be outside the fluid region is actually surrounded by the fluid. This figure confirms the physics behind bacterial propulsion in which a helix caused to rotate in a very low Reynolds number will necessarily translate (44
). A straight cylindrical tube (not shown) with the same physical conditions (e.g., viscosity of the surrounding fluid, stiffness constant, etc.) as in the helical filament except its shape yields no propulsion, as expected. (An animation movie that shows the propulsion of the model filament at a rotational frequency of f = 0.17
1 has been included in Supplementary Material.)
The total hydrodynamic force F on a flagellum can in general be obtained by integrating the force, f, per unit length or unit area, depending on the representation of the flagellum as either a line filament or a tube, respectively. Assuming ds to represent an element of length or area, in general, F can be written as
. The model described here represents the flagellum as a hollow tube. Thus, ds corresponds to an element of area. For a discrete representation, F can be obtained by summing the forces over the material points of the model. According to the SPH equation of motion in Eq. 5, the hydrodynamic force on each material point is the sum of the pressure force and the viscous forces arising from the interaction between the filament particles and the smooth fluid particles. The component of these forces along the helical axis adds up to provide a net thrust force that propels the model filament. The other components are expected to cancel each other because of symmetry. In Fig. 2, we show the x, y, z components of the hydrodynamic force for a frequency f = 0.017
1. The axis of the helix is chosen as the z axis, and the directions perpendicular to the z axis are assigned arbitrarily as the x and y axes. In the plots, the magnitudes of the instantaneous forces are shown together with the time-averaged values. The average forces in the directions perpendicular to the helix axis are found to be nearly zero, as expected, whereas a significant amount of force is observed along the long helical axis. This force corresponds to the thrust force. The evaluation of the average hydrodynamic force for different motor rotations (cf. Fig. 2 d) reveals that the magnitude of the thrust force increases with increasing
, whereas the magnitude of the other force components are still small.

View larger version (16K):
[in this window]
[in a new window]
|
FIGURE 2 Components of the total hydrodynamic force (a) along the filament axis (thrust force)(b) along an axis perpendicular to the helix axis, the x axis, and (c) along the y axis. The horizontal dashed line is a line marking the time average of each component. This is compared to a horizontal solid line marking the zero position. These lines are indistinguishable from each other for Fx and Fy. The forces are calculated for motor rotational frequency f = 0.017 1. (d) The magnitude of the time average hydrodynamic force components at different frequencies. The unit for the forces is given in reduced units.
|
|
Flow fields
One of the main advantages of this model, compared to, e.g., those using resistive force theory, is that the hydrodynamics is treated explicitly. This allows a more detailed examination of the behavior of flow fields around the flagellum. The recent PIV experiment on a macroscale model of bacterial flagella reveals the velocity distribution around rotating rigid helices and flexible helices (36
). For lack of similar numerical simulations done on flexible helices, the velocity fields from both the rigid and flexible helices are compared to a rigid helix based on the slender-body theory. The analysis in this section provides information on the nature of flow fields based on flexible helices, where a study is done on the flow fields generated by swimming filament. Note that the above studies are made using helices that are restrained from translating in contrast to a swimming filament of this model.
The flow around the model filament is studied by considering fluid particles over a plane bisecting the helical axis. First, the instantaneous velocities and pressures at each particle position rj are calculated. Then, the corresponding field variables at any space point r are measured by averaging the instantaneous values of all particles within a cutoff radius h of the space point. The averaging is done according to the SPH method, in which the contribution from each particle is weighted by the value of the kernel function at the particle's position. In solving the hydrodynamics, the boundary conditions need to be specified. Ideally, one wishes to surround the filament by an amount of fluid that is large enough to be regarded as an infinite fluid. However, this is limited by computational efficiency. In the model we describe, all external boundaries are fixed except the boundaries bisecting the filament axis. Periodic boundary conditions are applied along the filament axis. This allows the model filament to be surrounded by fluid at all times during propulsion, as mentioned above. The choice for a fixed external boundary on the remaining sides permits a comparison of the flow with known flows from fluid mechanics. Moreover, with this choice, our model attains qualitative similarity to the PIV experiment that is carried out inside a rectangular tank. It should be noted that other boundary conditions have also been tested, e.g., free surface sides. As expected, the flow pattern for the free surface is different from the fixed boundary. In the case of the free surface, the whole fluid rotates as a rigid body when the flow is fully developed. In the following, all analyses are based on the fixed boundary.
Velocity field
The velocity field at any point r is evaluated by SPH averaging of the velocities vj at positions rj that are within a cutoff radius of the point r, i.e.,
 | (9) |
The evaluation of the velocity field in this model is simple and straightforward, as compared to, say, the slender-body theory. In this method, this calculation amounts to summing individual particle velocity vj, weighted by the weighting function W, where the particle velocities are simply a result of integration of the SPH equation of motion. In contrast, the velocity field in the slender-body theory is found based on a linear superposition of the Stokeslets (point forces) and doublets (source dipoles) distributed along the centerline of the flagellum (49
). This involves the calculation of the Stokeslet and doublet tensors, as well as the strength of the point forces and source dipoles (19
,49
).
In Fig. 3, the calculation of the velocity vectors is shown for a fluid element over a cross section halfway between L/2 and L/2 in the z-direction, where L is the length of the simulation box. The analysis is done for flows generated by rotations of the model filament at motor speeds
= 2
x 0.0068
1 and 2
x 0.17
1. The values of the speed of sound corresponding to these frequencies are chosen such that M = 0.1. The viscosity parameter is set to
= 10 for both rotations. The Reynolds number for these selections is Re
0.006. The velocity vectors are shown for the times t = 3
and t = 5
for the rotational speeds
= 2
x 0.17
1 and
= 2
x 0.0068
1, respectively. These times correspond to the times when the flows are fully developed. In each of the figures, it is observed that a smooth flow of the fluid accompanies the rotation of the filament. This is expected for low-Re flows that are known to show a laminar behavior (50
). Notice that the velocity vectors are smoother for the lower speed. As expected, the laminar behavior becomes more evident for low speeds in the zero-Re limit. For other planes, patterns similar to those shown in Fig. 3 are observed. It is found that, unlike steady flows, the velocity vector deviates from a smooth flow pattern after long times; the point in time where this occurs depends on the frequency of motor rotation. The instability could be either due to the flexibility or the forward propulsion of the helical filament, or both. The flexibility may simply cause local disturbances that eventually propagate to other regions. The swimming of the model filament, on the other hand, may induce a shear stress that drags the fluid element in a direction parallel to the motion of the filament.
The magnitude of the instantaneous velocity field can be shown using 2D and 3D contour plots. Fig. 4 shows the contour plots of the velocity field for the rotation induced by a motor speed of
= 2
x 0.17
1. The 2D contour plot reveals the rotationally symmetric nature of the flow (except near the boundary). The 3D contour plot clearly demonstrates the decrease in the magnitude of the field for distances further away from the filament. These observations are typical of low-Re flow. To make a quantitative comparison between the flow induced by the model filament with a flow known from theoretical calculations, the average azimuthal velocity v
(r) is evaluated for the flow generated by the helix. The flow due to the model filament is an example of a flow with a moving surface. Noting the circular symmetry of the flow (neglecting the effect of boundary), the flow generated by the filament is similar to the Couette flow in the annular gap between two concentric cylinders of radii R and
R, one of them rotating at a rotational frequency
. For a steady flow, the velocity field of the Couette flow in cylindrical coordinates has the form v
= v
(r) and vr = vz = 0, where v
, vr, and vz are the azimuthal, r and z, components of the velocity, respectively. The analytical solution for v
(r) at any distance r from the center of the inner cylinder is given by (51
)
 | (10) |
To compare the theoretical value of v
with the results, the theoretical value from the above equation is fitted to the data. For this system, v
(r) is determined by radially averaging the velocity field evaluated by Eq. 9 over circles of constant radius r around the filament. The calculations are done for flow fields generated by different rotational frequencies
at the times when the flows are fully developed for each rotation.
Although the two systems are different in many respects, good qualitative agreement between the theory and the simulation is observed for the regions away from the filament. This qualitative similarity becomes more apparent when the theoretical value is plotted with the data, as shown in Fig. 5, where, for comparisons, the theoretical values are scaled down by a constant value of
0.35. Near the filament, a deviation is found from the qualitative behavior of the theoretical value. The origin of this deviation may appear to arise from the flexibility of the model filament studied. However, qualitatively similar behavior has been observed from the PIV experiment for both flexible and rigid rods (36
). In the PIV experiment, the flow field shows a sharp bend in the vicinity of the helix ring. A numerical simulation done on a rigid helix based on the slender-body theory showed the same behavior.
To understand the origin of this deviation for the present system, the azimuthal velocity is recalculated using the individual particle velocities, instead of the flow field derived from the SPH averaging (Eq. 9). The sharp decrease of v
(r) in the present system, as compared to the theoretical value, may arise from the velocity interpolation of Eq. 9. A better agreement between the theory and the simulation is found when v
(r) is fitted using the individual particle velocity (cf. Fig. 6). This leads us to conclude that, despite its flexibility, the flow due to the model filament appears momentarily as the Couette flow under a fixed external boundary.
Pressure field
The pressure field can be calculated in a manner similar to that for the velocity field. At each particle position ri, the pressure pi is evaluated using Eq. 8, based on the variations in density of the smooth particles. The pressure field p(r) at any point r in space is calculated by SPH averaging of the pressures pi as
 | (11) |
In Fig. 7 a, we show a 3D contour plot of the instantaneous pressure field for the flow shown by the velocity field in Fig. 4. It is interesting to note that the contour plot of the pressure field looks like a flipped velocity field. That means the pressure field is maximum where the velocity field is minimum. Note that this is different from the pressure field of a Stokes flow in an infinite fluid. In the Stokes steady flow, the pressure field of a singular point force concentrated at r = 0 has a radial component that decays like 1/r2, yielding a vanishing pressure at infinity. Since the Stokes equation is linear, the superposition of such singular points will also behave similarly. This difference is a result of the boundary condition used.
As in the case of the velocity field, the behavior of the pressure field depends on the boundary condition set to the fluid. For the fixed external boundary, following an argument similar to that leading to Eq. 10, an expression for the pressure gradient along r can be found as
 | (12) |
This is nothing but the centrifugal force acting on a fluid element of mass density
. The pressure gradient provides the force needed to maintain the flow. Consequently, the velocity is large in a region where
p/
r is large, and small otherwise. Theoretically speaking, if the above equation is integrated using the expression for the azimuthal velocity (Eq. 10), it can be observed that p(r) is maximum at the external boundary. Fig. 7 b confirms this argument. In the figure, we plot the radially averaged pressure field obtained from the data of Fig. 7 a. The figure also shows the gradient of the average pressure field. It is evident from the figure that the speed of the flow is large where the pressure gradient is large.
Another interesting feature of the pressure field is the early-stage development of the field, i.e., during the time before the flow is fully developed. Fig. 8 shows a surface plot of the pressure field evaluated at different times, for the fluid element studied in Fig.7. The figure depicts the time evolution of the pressure field as the flow stabilizes. When the model filament begins to rotate, the helix induces a pressure gradient caused by an increase in the density of the closest layer. This creates a pressure build-up around the filament. With time, the nearest layer influences its surrounding due to the pressure gradient. This effect spreads out to other layers, causing the pressure to affect regions further from the filament until the whole system stabilizes, in a manner similar to a wave propagating in a medium. It should be noted that in a real bacterial motion one observes a series of stops and runs, where the cell stops momentarily during the tumbling process. Thus, the early-stage development seen here is not an isolated event that only happens at the beginning of the entire chemotaxis. It is rather part of the bacterial propulsion process, although in a more complex manner in the real system.
Deformation: rotational instability
During the tumbling process, when one or more rotary motors reverse their directions, the flagella undergo deformations that involve a complex sequence of shape transformations. When such polymorphic transitions occur, a range of motor speeds is observed that may differ from the speeds when the cell runs smoothly. Numerical studies have been performed on elastic filaments to study the nature of deformations arising from varying rotational speeds (32
,52
). These studies reveal two dynamical regimes of motion depending on the rotation rate. They find that when an elastic rod is rotated at a particular rotational speed
, the shape of the rod becomes unstable depending on the speed of the rotation. At low rotational speeds, the rod spins about its long axis, with the centerline remaining straight. This is referred to as twirling. When the speed is high, the rod rotates in a form where the centerline is bent, referred to as whirling (52
) or overwhirling (32
). These studies are done on straight rods in a viscous hydrodynamic environment. More complex phenomena are expected with a rotating helical filament. In this section, the dynamical instability for the model helical filament is examined. The focus is on the influence of hydrodynamics on the deformation behavior of model flagella. To clearly identify the effect of viscous interactions, the nature of deformation of the model filament is compared with and without hydrodynamics. The model discussed here is capable of this extension, in contrast to, for example, the immersed-boundary method (32
). The immersed-boundary method is inherently designed to incorporate nonlocal hydrodynamic effects. On the other hand, in Wolgemuth (52
), hydrodynamics is treated at the local level in terms of drag forces.
Consider an isolated helical filament (without hydrodynamics) that is rotated at one end with the other end free, under a boundary condition where the rotated end is pinned or restrained from translating. As the base of the filament is caused to rotate at a particular frequency
, the entire structure responds to the torque in accordance with the stiffness constant and the speed of motor rotation. As in other studies, a twirling to whirling transformation is observed with increasing
. In Fig. 9, the time sequence of the helix rotation is shown for the rotational frequencies f = 0.017
1, 0.17
1, and 0.85
1. At f = 0.017
1, the centerline of the filament remains straight at all times, indicating that twirling is stable. For f = 0.17
1 and 0.85
1, twirling is unstable and the centerline of the filament bends upon rotation, with the deflection being large for the higher frequency. This observation is similar to the whirling of shafts described in Love (53
). Furthermore, the shape of the instability resembles those described as crankshaft and speedometer motion in the analytical study of bistable helices (25
).

View larger version (79K):
[in this window]
[in a new window]
|
FIGURE 9 Twirling and whirling motion of an isolated helical filament pinned at the base and rotated at frequencies (a) f = 0.017 1, (b) f = 0.17 1, and (c) f = 0.85 1. The different colors represent the same filament at different times. For a, the times are t = 1 (red), t = 50 (cyan), and t = 92 (blue). For b, the times are t = 1 (red), t = 6 (cyan), and t = 10 (blue). For c, the times are t = 1 (red), t = 2.5 (cyan), and t = 5 (blue).
|
|
The effect of hydrodynamics may be understood by repeating the above calculation for a filament surrounded by a fluid with a viscosity parameter
= 10 and speed of sound adjusted to yield a low Mach number M = 0.1. As in the isolated helical filament, no deflection is observed for f = 0.017
1. However, it is also observed that no deflection occurs for f = 0.17
1. This is in contrast to the case of an isolated model filament, where a small deflection is found for this speed of rotation. The lack of instability for f = 0.17
1 indicates that hydrodynamics plays a role in suppressing small deflections. For sufficiently high rotational speed, such instabilities become inevitable. In this case, the nature of the deformations depends on another factor as well, i.e., the deflection depends not only on the rotational speed but also on the viscosity of the fluid. It should be recalled that the SPH model does not provide a direct measure of viscosity since it is designed based on an artificial viscosity. However, from the expression for Re that can be derived on the basis of scaling analysis, it becomes clear that
can be related to
and c by
= 
ch. For a given speed of sound, the viscosity parameter
can be varied to study qualitatively the effect of viscosity.
Fig. 10 shows the deflections at high rotational speed of f = 0.85
1 for the values of
ranging from 0.01 to 10. For this calculation, the criterion for M is relaxed slightly to achieve stable simulation. For the motor rotational frequency of f = 0.85
1, the speed of sound is chosen as c = 55 (in reduced units) such that M changes from M = 0.1 to M = 0.17. This change results in a <3% density variation. The overall qualitative behavior shown in the plot reveals the progressive change of the deformation behavior with increasing viscosity. For
= 0.01, a deflection similar to the isolated helical filament is observed, where the centerline of the filament bends while rotating. This indicates that twirling is unstable for this viscosity. With increasing
, the deflection decreases progressively, eventually vanishing at
= 1. For high viscosity,
= 10, the deflection is replaced by a deformation that winds part of the filament around itself. This deformation shows the effect of rotary stress in changing the helical pitch or radius of the filament. Intuitively, this is what is expected from an elastic object twisted in a highly viscous medium. On the other hand, the observation hints at the effect of mechanical stress in leading to polymorphism. Note that a polymorphic transition in a real flagellum is a complex process that involves chirality reversal, where a sequence of transformations from normal to supercoiled to curly-1 and then to normal is typically seen during a single run-and-tumble cycle (38
). This model captures only the factor that may trigger this complex process. Future improvements are needed to capture the entire polymorphism through a model that fully incorporates both hydrodynamics and elasticity.
 |
CONCLUSIONS
|
|---|
A simple hybrid model based on the elastic network model and the smooth-particle hydrodynamics methods has been developed to study the propulsive dynamics of bacterial flagella in a viscous hydrodynamic environment. The propulsion of the bacterial flagellum involves the interplay of elasticity and hydrodynamics. The model incorporates these two aspects by explicitly modeling the elastic nature of the filament and the hydrodynamics of the surrounding fluid. The model further incorporates the structural organization of the flagellum at a coarse-grained level in terms of the 11 protofilaments based on the physical properties of flagellin, the protein building block of flagellum. The overall elastic material property of the real flagellum is matched to the stiffness constant of the model filament by parameterizing the model using the flexural rigidity of flagellum found from experiments. The elastic helical filament representing flagellum is coupled to the hydrodynamics through a boundary condition that enforces a no-slip boundary condition. The effect of the surrounding fluid on the model filament is manifested through the viscous and pressure forces. In the future, an explicit multiscale connection of this model to atomistic or coarse-grained molecular dynamics simulations may be possible.
The physics of bacterial propulsion has been captured by this simple model through the coupling of the helical filament and the surrounding fluid. The model also allows for a rigorous account of flows surrounding the filament. To study the hydrodynamic properties, a fixed boundary condition is set to those sides parallel to the long helical axis. The instantaneous velocity field resulting from the rotational motion of the model filament shows a behavior expected from low Re flows. The velocity vector shows a laminar flow and the magnitude of the velocity field decays as the distance from the filament. For a fully developed flow, the qualitative feature of the velocity field is similar to that observed from the PIV experiment and the slender-body calculation. The far-field behavior of the velocity field is qualitatively similar to the Couette flow, where a good fit to the data is found when the theoretical value is scaled down for comparison. The individual particle's velocity fits well with the theoretical predictions of Couette flow for all regions upon a similar scaling. The forward propulsion of the filament eventually causes a disturbance on the steady flow. The pressure field of a fully developed flow is shown to maintain the flow by providing a pressure gradient that is relatively high near the filament. The pressure gradient acts as a centrifugal force. The development of the pressure field before the flow is fully developed has also been investigated. The pressure disturbance propagates outward in a manner similar to a wave.
The deformation of the model filament was studied in the presence and absence of the surrounding fluid. In the absence of the surrounding fluid, the filament undergoes a dynamic instability where the rotation of the helix changes from a twirling state to a whirling state. The deflection of the filament increases with the increase of rotational speed. It is found that the hydrodynamic effect via the SPH solvent suppresses the deflection resulting from small motor rotations. For high motor rotations, the deformation depends on the viscosity of the fluid. For relatively low viscosity, the filament shows similar dynamic instability as that found for an isolated filament. Upon increasing the viscosity, part of the filament twists around itself so as to change the shape of the helix. The complex process of polymorphism that involves shape transformations and chirality reversal may be initiated due to mechanical loadings, as qualitatively observed in the present system. To capture the full polymorphic transformation, some modifications are needed to the model presented here. In particular, a model of a flagellum that incorporates some level of atomistic detail may be required. Another aspect of our model that requires future improvement corresponds to the artificial viscosity of the SPH model. This factor limits the capability of the model to do quantitative comparisons with other theoretical or experimental predictions. Future improvements might be to employ another variant of the SPH method known as smooth-particle applied mechanics (SPAM) (54
,55
). The latter method is the same as the SPH method except for the artificial viscosity.
 |
SUPPLEMENTARY MATERIAL
|
|---|
An online supplement to this article can be found by visiting BJ Online at http://www.biophysj.org.
 |
ACKNOWLEDGEMENTS
|
|---|
This work was supported by National Institutes of Health grant No. GM053148.
Submitted on June 14, 2006;
accepted for publication August 11, 2006.
 |
REFERENCES
|
|---|
1. Berg, H. C. 2004. E. Coli in Motion. Spring-Verlag, New York.2. Berg, H. C. 2000. Motile behavior of bacteria. Phys. Today. 53:2429.
3. Namba, K., and F. Vonderviszt. 1997. Molecular architecture of bacterial flagellum. Q. Rev. Biophys. 30:165.[CrossRef][Medline]
4. Berg, H. C., and R. A. Anderson. 1973. Bacteria swim by rotating their flagellar filament. Nature. 245:380382.[CrossRef][Medline]
5. Silverman, M., and M. Simon. 1974. Flagellar rotation and the mechanism of bacterial motility. Nature. 249:7374.[CrossRef][Medline]
6. Asakura, S. 1970. Polymerization of flagellin and polymorphism of flagella. Adv. Biophys. 1:99155.[Medline]
7. Calladine, C. R. 1983. Construction and operation of bacterial flagella. Sci. Prog. 68:365385.[Medline]
8. Calladine, C. R. 1975. Construction of bacterial flagella. Nature. 255:121124.[CrossRef][Medline]
9. Samatey, F. A., K. Imada, S. Nagashima, F. Vonderviszt, T. Kumasaka, M. Yamamoto, and K. Namba. 2001. Structure of the bacterial flagellar protofilament and implications for a switch for supercoiling. Nature. 410:331337.[CrossRef][Medline]
10. Hasegawa, K. I., I. Yamashita, and K. Namba. 1998. Quasi- and nonequivalence in the structure of bacterial flagellar filament. Biophys. J. 74:569575.[Abstract/Free Full Text]