Chirality-induced bacterial rheotaxis in bulk shear flows

Mean rheotactic velocity

A dilute suspension of E. coli bacteria is injected at a given flow rate Q into a microchannel of width W = 600 μm and height H = 100 μm (Fig. 1A), which imposes to a very good approximation a planar Poiseuille flow V_x(z) = 4V_maxz(H − z)/H² sufficiently away from the side lateral walls, with V_max, the flow velocity in the center of the channel. The local shear rate is then γ.=∂Vx(z)/∂z=4Vmax(H−2z)/H2. Using a high-magnification objective lens, bacteria tracks are recorded as projections into the x-y plane at different distances from the bottom wall (Fig. 1B) and far away (≳250 μm) from lateral side walls. Exact calibration of the local flow velocities and channel orientation is obtained by the simultaneous recording of passive tracers for all experiments performed. This calibration step is essential for the determination of bacteria orientation and velocity distributions. Since we focus on bulk dynamics, bacterial trajectories close to the top and bottom channel walls are not included in our analysis.

• Download high-res image
• Open in new tab
• Download Powerpoint

Figure 1 (C and D) show examples of bacterial trajectories in layers in the upper and lower parts of the microchannel Poiseuille flow. One clearly sees a bias toward the right with respect to the negative flow direction on (C) and toward the left on (D), confirming the rheotactic cross-streamline migration in the y direction induced by the chirality of the left-handed bacteria flagella (37) and being a function of the sign of the local shear rate.

Quantitative measurements of mean bacteria velocities (v¯x,v¯y) and transport velocities of passive tracers (V_x,V_y) for different flow rates Q are displayed in Fig. 2 (A and C) in the flow direction (x direction) and in the vorticity direction (y direction), respectively, as a function of the distance z from the bottom wall. The passive tracers (open symbols) follow a Poiseuille profile [V_x(z) in Fig. 2A] without a drift in the y direction (Fig. 2C), indicating the perfect alignment of channel and microscope, and the x axis is identical to the flow direction.

• Download high-res image
• Open in new tab
• Download Powerpoint

Different local shear rates γ. can be obtained in two ways, either by varying V_max via the imposed flow rate Q or by varying the position z inside the channel. For a systematic analysis of bacterial rheotaxis, allowing to independently investigate the role of local shear rate and positions within the Poiseuille flow, both Q scans and z scans are performed.

Mean bacteria velocities (filled symbols) deviate from the background velocities both in the parallel (x) direction and in the transverse (y) direction. In the x direction, bacteria velocities are higher in the center of the channel and lower close to channel boundaries compared to the background flow. This is illustrated more clearly in Fig. 2B where the z-dependent relative velocities between bacteria and background flow are shown. This indicates a preferred orientation of bacteria downstream in the channel center and upstream closer to the channel walls for all considered flow rates. Mean bacteria velocities v¯y confirm the visual observations from Fig. 1 (C and D) and have opposite signs in the lower and upper half of the Poiseuille flow. This rheotactic velocity v¯y increases with increasing flow rate Q and with decreasing distance z from the wall. Note that wall effects are visible for the data points at a distance of 0.1H (∼10 μm) from both top and bottom walls and will be excluded from further analysis in this paper.

Figure 2C shows mean drift velocities v¯y as a function of local shear rates from different datasets, including the results from the three z scans from Fig. 2B and from two different Q scans. We note that we scale the velocities by the average bacteria velocity v₀ = 25 μm s⁻¹ in the absence of flow (see the Supplementary Materials). Representing the mean velocities as a function of the local shear rate leads to a reasonable data collapse, indicating that the local shear rate is the main control parameter of our system. Within our range of shear rates, the increase in mean drift velocity is first linear with local shear rate, similar to what was observed in (37), and then reaches a maximum of around half the average bacteria swimming velocity at the highest shear rate.

Theoretical framework

To understand the physical mechanisms behind bacterial rheotaxis in microchannel flow, we develop a theoretical framework that captures the dynamics of individual noninteracting bacteria. We account for the elongated and chiral shape of the flagellated bacteria, their self-propulsion velocity v, tumbling, translational and rotational noise, and their advection and rotation in Poiseuille flow. Denoting the instantaneous position of a bacterium by r and its orientation by e, we can write the dynamics of a bacterium in Poiseuille flow asdrdt=ve+Vx(z)+ℋ·ξ(1)for the positional dynamics anddedt=[ΩJ+ΩC+2Drξr]×e(2)for the orientational dynamics. For Poiseuille flow, these two equations are coupled, and Eq. 1 includes orientation-dependent self-propulsion, advection in Poiseuille flow, and anisotropic translational diffusion, and Eq. 2 includes position-dependent reorientation in flow. The random numbers ξ_i and ξir model Gaussian white noise (see Materials and Methods); ℋ is related to the anisotropic translational diffusion tensor (see Materials and Methods), and D_r is the rotational diffusion constant. We use the well-known Jeffery reorientation rate Ω^J to capture the rotation in flow due to elongation (9, 20, 21, 39)ΩxJ=γ.12GexeyΩyJ=γ.12(1+G(ez2−ex2))ΩzJ=−γ.12Geyez(3)with the Bretherton shape factor G = (α² − 1)/(α² + 1) ≲ 1, where α is the effective aspect ratio of the bacteria. The chiral strength ν of a bacterium is related to the specific shape of the left-handed chiral flagella bundle and the size of the head (37, 38), leading to a chirality-induced reorientation rate Ω^C, as derived in (38)ΩxC=−νγ·ez(2ex2−1)ΩyC=−2νγ·exeyezΩzC=−νγ·ex(2ez2−1)(4)

The same result has been rederived recently in (40). To model tumbling, we include tumble events that modify the bacterium orientation instantaneously at exponentially distributed times (see Materials and Methods). We note that we neglect a term in Eq. 1 that captures a small passive rheotactic drift force but is negligible compared to the other terms for our regime of parameters.

The theoretical model (Eqs. 1 to 4) contains a number of parameters that need to be determined. The swimming speed v of individual bacteria is assumed to be normally distributed with mean v₀ = 25 μm s⁻¹ and SD of 8 μm s⁻¹ to closely match experimental conditions in the absence of flow (see the Supplementary Materials). While the value of D_r = 0.057 s⁻¹ (41) and the tumbling statistics (42) can be estimated from the literature, we treat α and ν as free parameters. These parameters can be adjusted using the precise experimental results. We thus perform Brownian dynamics simulations of the coupled Eqs. 1 and 2 including tumbling and simple steric repulsion for swimmer-wall interactions to be compared to the experimental results for the mean rheotactic velocities v¯y as a function of the local shear rate γ. as shown in Fig. 2C.

The initial slope of v¯y(γ.) can be used to adjust the mean chiral strength of the bacteria to ν = 0.06, being of the same order as estimated previously using resistive force theory (37, 38). The parameter ν only depends on the bacterial shape. In particular, its value depends nontrivially on the helical shape of the flagella bundle, is maximized for large cell bodies, and vanishes in the limit of very small cell body size. The effective aspect ratio α is determined from the deviation from the linear regime toward the saturation of the rheotactic velocity at high shear rates, with the results not being very sensitive to the value of the parameter α. Here, we use α = 5, in agreement with typical experimental values (21). These values for α and ν are used for all further comparison between experiments and simulations.

From Fig. 2D, we can see that by adjusting these parameters, the numerical model describes very well the experimental observations. As in the experiments, data obtained from different positions z in the channel and different flow strength collapse onto almost a single curve, but small deviations stem from the velocities obtained from different layers in the channel.

In Fig. 2B, we compare the longitudinal velocity differences between bacteria and the Poiseuille flow from numerical simulations with the experiments, which show the same results, i.e., a preferred upstream orientation near the walls and a downstream orientation in the center. This is in excellent agreement with a kinematic model, which neglects the swimmer chirality (22). This indicates that chirality does not influence, at least qualitatively, the polarization profile along the flow direction but mainly affects the bacteria orientations in the transverse direction.

Velocity and orientation distributions

We now turn to velocity and orientation distributions. Figure 3A shows the rheotactic velocity distributions P(v_y) at a given height z = H/4 for different local shear rates corresponding to the mean values shown on Fig. 2C. As expected, without flow and for small shear rates, the velocity distribution is very broad and centered around zero. For increasing shear rates, the distributions are shifted more and more toward positive v_y values, together with a small decrease of the width of the distribution. The experimental and numerical distributions are in excellent agreement and reinforce the validity of our model. Note that using a swimming velocity distribution matching the experimental results is crucial to obtain such an agreement, whereas it is sufficient to work with average values for α and ν.

• Download high-res image
• Open in new tab
• Download Powerpoint

The rheotactic velocity of each bacterium is a direct consequence of its three-dimensional (3D) orientation and its intrinsic swimming velocity v setting the instantaneous velocity v_y(t) = ve_y = v cos Θ(t) sin Ψ(t). Here, the angle Θ is linked to the orientation of the bacterium in the z direction, e_z = sin Θ, and e_x = − cos Θcos Ψ, and e_y = cos Θsin Ψ (see Fig. 1B). While bacteria orientations can directly be determined from the numerical simulations, only projections of bacteria trajectories into the x-y plane are accessible from the experiments and the in-plane orientations are obtained from the orientation of the velocity vector (see Materials and Methods), defined by the angle Ψ. We thus show the corresponding orientation distributions P(Ψ) for different local shear rates in Fig. 3B for the layer z = H/4. Figure 3C shows a sketch of the bacterium orientation Ψ and relation to up- and downstream and left and right orientation. At zero shear rate, the orientation distribution is flat and remains so for small shear rates, indicating no preferred bacteria alignment under weak flow. With increasing shear rate, bacteria orient, on average, more and more toward positive values of Ψ and thus toward the right with respect to the negative flow direction. In addition, a double peak is emerging, which is found to be mostly symmetric around Ψ = π/2, corresponding to a bacteria orientation perpendicular to the flow direction. Bacteria are thus oriented up and downstream around the perpendicular orientation. These peaks become more pronounced and also move closer together with increasing shear rate. They correspond to a single peak in the velocity distributions as seen on Fig. 3A.

To represent a more complete picture of the orientation distributions under flow, Fig. 4 shows the distributions layer by layer along the channel height. z scans have been performed at different imposed flow rates Q and are represented as a function of the wall shear rate γ.W=4Vmax/H. Higher flow rates are displayed only for the simulation results, as at those large flow velocities, the time resolution of the experimental image capture is not sufficient to resolve the bacteria trajectories. The probabilities of the rheotactic velocity and orientation distributions are represented with a color code where yellow corresponds to a high probability and blue corresponds to low probability. Peaks in the distributions are thus easily identified as regions of yellow color. For low flow rates, the continuous shift of the peak velocity from zero toward higher velocities when moving away from the middle of the channel toward the channel walls (and thus with increasing local shear rate) is clearly visible, in agreement with Fig. 2C. For larger flow rates, in layers closer to the channel walls and thus for higher local shear rates, only a very weak increase or even a saturation of this peak value can be identified, in particular, from the numerical results. From the orientation distributions, the existence of a double peak is clearly visible for moderate flow rates and a decrease of the distance between the two peaks is observed when getting closer to the channel walls. At high flow rates, the simulation results show that the double peak merges into a single but wide peak.

• Download high-res image
• Open in new tab
• Download Powerpoint

These figures also show the up- and downstream orientations induced by the swimming dynamics of the bacteria in the Poiseuille flow, where bacteria cross multiple layers in the z direction during a trajectory (19, 21, 22, 43). This is, in particular, visible at small wall shear rates, where in the center of the channel, most bacteria are oriented downstream, as can be seen by the yellow peak close to π. Closer to the channel walls, this peak is found around 0, corresponding to an upstream orientation, in agreement with the relative bacteria velocities shown in Fig. 2B.

Theoretical understanding of the observed orientation distributions

To explain the orientation distributions P(Ψ, Θ) at the origin of the rheotactic behavior, we start by discussing a simple model system for bulk rheotaxis, namely nontumbling elongated bacteria in linear shear flow with constant shear rate γ.. In the absence of noise, the equations for the orientations of chiral microswimmers (Eq. 2) can be rewritten in terms of the angles Ψ and Θ and read [see also (38)]Ψ.=γ.2(1+G)sin Ψ tan Θ+γ.ν cos Ψcos 2Θcos ΘΘ.=γ.2(1−G cos 2 Θ)cos Ψ+γ.ν sin Ψ sin Θ(5)

Both the Jeffery (first terms) and the rheotactic contributions (second terms) are linear in the shear rate γ., which, as a consequence, only determines the time scale of the dynamics.

The orientation (Ψ(t), Θ(t)) of nonchiral swimmers (ν = 0) simply follow Jeffery orbits, which are solutions of Jeffery’s equations of motion (see also Eq. 3). The Jeffery orientation phase space is shown in Fig. 5A, and closed streamlines indicate the well-known periodic Jeffery solutions that depend on the initial orientation of the swimmer and correspond to a constant of motion, the so-called Jeffery constant (9). The arrows indicate the direction of motion on the unit sphere. Again, the shear rate γ. does not alter the Jeffery orbits but only sets the time scale of how fast an elongated particle or bacterium rotates. Note that when a particle starts oriented to the left/right (Ψ < 0/Ψ > 0), it never switches to the other side.

• Download high-res image
• Open in new tab
• Download Powerpoint

Adding the rheotactic term now breaks the left-right symmetry, and as can be seen in Fig. 5B, all orientations eventually end up on the right (Ψ > 0), leading to the mean rheotactic velocity v¯y>0. To understand the orientational behavior in more detail, we see that (Ψ* = ± π/2, Θ* = 0) are two fixed points of the system (Eq. 5), defined where (Ψ.*,Θ.*)=(0,0). The eigenvalues of the stability matrix for the linearized system around the fixed points can be evaluated toλ1,2=±i2γ.1−G2−4ν2=±iω(6)

These Eigenvalues have a vanishing real part, which means that the fixed points are actually not stable/unstable but marginally stable. In other words, the typical time scale at which a bacterium orients toward the fixed point on the right diverges, and it will therefore never be reached. We observe that simulated trajectories without noise do not reach the fixed point but get trapped in periodic orbits around the marginally stable fixed point.

We now turn to the question of how rotational noise, always present in experiments, influences the orientational distribution of bacteria in flow. The rotational diffusion constant D_r of the bacterium compared to the strength of the shear rate γ. now plays an important role. For nonchiral swimmers (ν = 0), noise affects the dynamics similar to passive elongated particles (44). Because of noise, trajectories now erratically move around between different Jeffery orbits (12), which strongly influences the orientation distribution function P(Ψ, Θ), as shown in Fig. 5C. Here, we show results of 1000 averaged steady-state trajectories started uniformly distributed on the unit sphere for a fixed rotational diffusion D_r = 0.057s⁻¹ but for a broad range of different shear rates. Purple and blue correspond to low probability, while yellow and red correspond to high probability (see scale bar). The first row in Fig. 5C shows the case where noise is not yet important because of the high shear rate γ.=1000s−1, and P(Ψ, Θ) almost follows the deterministic solution obtained from Eq. 5. The peaks at (Ψ = 0, Θ = 0) and (Ψ = π, Θ = 0) correspond to particles aligned with the flow, upstream or downstream, respectively. These peaks are a consequence of the fact that elongated particles moving on a (typically kayaking) Jeffery orbit have a high probability to come close to these positions, which can be seen by the increased density of streamlines in Fig. 5A. Together with a decreased rotation rate near these points, this leads to the large probability observed. Note that the slowest rotation is near the log-rolling states (Ψ = ± π/2), but hardly any initial condition leads to orbits that come close to these states and precludes high probabilities around them.

As can be seen in the second row in Fig. 5C, reducing the shear rate by a factor of 10 does not change the distributions substantially, meaning that noise is still not important for nonchiral swimmers. This changes for smaller shear rates, γ.=10 s−1, where the peaks of the distributions become less pronounced. At even smaller flow rates, γ.=1 s−1 (last row), noise randomizes the orientations almost completely, and only very weak peaks remain. Note that both upstream and downstream peaks are not symmetrically distributed around Θ = 0 for small γ. because of a subtle interplay between noise and flow, which leads to nonuniform microswimmer profiles in the z direction (22, 43, 44). In general, the important physical quantity that determines the nonchiral orientation distribution is the rotational Péclet number Per=γ./Dr. As we will see in the following, the situation is less simple when particles are chiral.

For chiral microswimmers (ν > 0), the signatures of the marginally stable fixed point can be seen at high shear rates. As shown in the first row of Fig. 5D, even for very weak noise or high shear rates (γ.=1000s−1) the peak of the distribution is not a simple δ peak but smeared out around the marginally stable fixed point at Ψ = π/2. For smaller shear rate γ., when noise becomes more dominant, the probability distribution becomes even broader and develops bimodal peaks, which shift more and more toward orientations aligned (or antialigned) with the flow when lowering γ. further (second and third row in Fig. 5D). The positions of these peaks are triggered by a competition between the chirality-induced attraction to the right and the Jeffery peaks, which become more and more dominant for smaller shear rates, since noise then helps to move around in the phase space more easily. For low γ. (last row), chiral migration to the right is very inefficient compared to noise, and the probability distributions look very similar to the nonchiral system.

Last, we note that the frequency of oscillation ω (see Eq. 6) is comparable to the pure Jeffery frequency ωJ=γ.1−G2/2 but is reduced by approximately 5% for our standard parameters ν = 0.06 and α = 5.

Universal scaling and master curve

We now analyze the dependence of the rheotactic drift on the system parameters. Therefore, we simulate the dynamics of different rotational diffusion D_r, bacterium aspect ratio α, and chiral strength ν for varying shear rates γ.. Figure 6A shows the scaled mean rheotactic velocities v¯y/v0 as a function of the shear rate for different nontumbling bacteria in simple shear flow. As expected, none of the curves approach v¯y/v0→1 even for very large shear rates because of the aforementioned trapping and oscillations of the swimmer orientations in the vicinity of the marginally stable fixed point (Ψ* = π/2, Θ* = 0).

• Download high-res image
• Open in new tab
• Download Powerpoint

As seen before, for small shear rates, all curves increase linearly with the shear rate (Fig. 2D). The reason is that for small shear rates, the peaks of the distributions are close to Ψ^± = {0, π} and Θ^± = 0 (see also Fig. 3B), and the peaks linearly shift away from them with increasing shear rate. The slope is stronger for larger chirality ν or weaker rotational Brownian noise D_r.

To understand the behavior at small shear rates, we linearize Eq. 5 around (Ψ^±, Θ^±), i.e., evaluate it at angles Ψ = Ψ^± + ϵ_Ψ and Θ = Θ^± + ϵ_Θ with small ϵ_Ψ and ϵ_Θ. To the lowest order in ϵ_Ψ and ϵ_Θ, we then obtain Ψ.(Ψ,Θ)=±γ.ν + h.o.t. and Θ.(Ψ,Θ)=ΩΘJ + h.o.t., where all higher-order terms (“h.o.t.”) are at least quadratic in the angles ϵ_Ψ and ϵ_Θ and ΩΘJ is the Jeffery contribution to the angular velocity component Ω_Θ. The reorientation rate Ψ., which is responsible for the drift toward the right, is to the lowest order independent of G and linear both in γ. and in ν. This suggests that the shifts of the peaks Ψ* of the probability distributions, and, hence, of the mean values v¯y, should be a function of γ.ν, acting against the noise quantified by D_r.

Our analysis and the results shown in Fig. 6A now suggest introducing a dimensionless chirality numberC=γ·νDreff(7)as the main relevant physical quantity, which regulates rheotaxis. Here, we use Dreff as the effective rotational diffusion constant determined via the exponential decay of the bacterium orientation autocorrelation function in the absence of flow, 〈e(0)·e(t)〉=e−2Drefft (23). It is simply the rotational diffusivity D_r for nontumbling bacteria but larger for tumbling bacteria because of the extra tumbling reorientations and estimated to Dreff≈0.25 s−1 in our experiments. Note that together with the experimentally explored range of local shear rates γ.≈0 to 50 s−1, this corresponds to a range of Péclet number Pe ≈ 0 to 200.

In Fig. 6B, we now show the same data as in Fig. 6A but as a function of the chirality number 𝒞. As can be seen in the inset of Fig. 6B, for small shear rates, all curves fall on top of each other. At higher shear rates, only curves with the same α (and hence G) fall on top of each other.

This result is important, as it points in this weak shear regime, on the relevant parameters defining the rheotactic drift velocity for a microswimmer. The linear relation between mean transverse rheotactic velocity v¯y and shear rate γ., v¯y=l0γ., also reveals a length l₀ quantifying how far a bacterium drifts in the vorticity direction at unit shear rate. For our bacterial strain, this length is l₀ ≈ 0.75 μm (see Fig. 2D), and we now know from Eq. 7 and the prefactors shown in the insets of Fig. 6 (B and C) that l₀ is related to the bacteria parameters via l0≈0.1v0ν/Dreff. We note that the prefactor of ∼0.1 is a purely numerical constant that does not depend on the system parameters.

At larger shear rates, the bacteria with lower aspect ratio (α = 3; diamonds) have a higher rheotactic velocity, and the swimmers with higher aspect ratio (α = 10; squares) have a lower rheotactic velocity. This can be understood by looking at higher orders of the stability around the locally marginally stable fixed point around Ψ = π/2. It can be shown that away from the linearly marginable stable fixed point, the attraction toward this point is stronger for smaller α and, hence, the rheotactic velocity is larger.

We also test our scaling for tumbling bacteria and compare the results to nontumbling bacteria in Fig. 6C. We can see that the scaling also works well for tumbling bacteria, where all of the data (black symbols) collapse onto a single curve. The nontumbling and tumbling master curve do not exactly fall on top of each other, which suggests that the simple scaling with Dreff is not perfect. Note that, as described before, we include tumbling not only as defining an enhanced effective rotational diffusion but also by adding random tumbling events on top of the continuous rotational diffusion. While, in the absence of flow, long-time dynamics may be captured by introducing an effective temperature from stemming from an effective diffusion constant (23), this concept fails to be the accurate measure in our situation, where the interplay of flow and tumbling prohibits a simple effective temperature scaling. Again, we obtain a linear regime for small γ. (inset) where the slope is almost the same (but not exactly) as that for the nontumbling bacteria.

Now, we turn to the case of tumbling bacteria in Poiseuille flow by plotting the data extracted from different layers in the channel (and hence local shear rates), as shown in blue in Fig. 6C. The fact that this data collapses together with the simple shear data on a single curve reiterates the fact that, at least for our channel height (H = 100 μm), bacterial rheotaxis is determined by the local shear rate.

Last, we show in Fig. 6D how the curves in Fig. 6C approach toward their asymptotic values. As discussed before, bacteria approach the infinite shear limit very slowly, which can be fitted, as a first approximation, to a power law ∼𝒞^−0.5, as a consequence of the marginally stable fixed point at Ψ = π/2.

Comparison to experimental and numerical results

Our theoretical analysis allows us now to understand all the observed aspects of the rheotactic velocities in both experiment and Brownian dynamics simulations. First, the aforementioned result that the position of the peaks in P(Ψ) for chiral bacteria move closer and closer together with increasing shear rate (Fig. 5D) is clearly observed in the experiments and simulations in Poiseuille flow using tumbling bacteria, as shown for P(Ψ) in Figs. 3B and 4B.

Second, we can now understand why the plateau value of the mean rheotactic velocities obtained in experiments and simulations for high shear rates does not approach the mean free swimming speed v₀ (Fig. 2C): The occurrence of a marginally stable fixed point traps bacteria in periodic-like orbits, resulting in broad orientational distributions.

Third, we are able to predict the linear increase of the mean rheotactic velocity for small shear rates, as seen experimentally and in simulations in Fig. 2C, and to determine the dependence of the mean rheotactic velocities on the relevant parameters, resulting in a universal scaling law.

Fourth, we demonstrate that our analysis in simple shear flow qualitatively explains the behavior in Poiseuille flow. We find that the orientation probability distributions P(Ψ, Θ) in Poiseuille flow obtained from layers with local shear rate γ. follow the respective distributions in simple shear flow. Tumbling allows bacteria to exploit nonpopular regions in phase space more frequently compared to nontumbling bacteria, but the qualitative results are comparable to the ones of nontumbling bacteria.