Fluidic bacterial diodes rectify magnetotactic cell motility in porous environments

Directed motility enables swimming microbes to navigate their environment for resources via chemo-, photo-, and magneto-taxis. However, directed motility competes with fluid flow in porous microbial habitats, affecting biofilm formation and disease transmission. Despite this broad importance, a microscopic understanding of how directed motility impacts the transport of microswimmers in flows through constricted pores remains unknown. Through microfluidic experiments, we show that individual magnetotactic bacteria directed upstream through pores display three distinct regimes, whereby cells swim upstream, become trapped within a pore, or are advected downstream. These transport regimes are reminiscent of the electrical conductivity of a diode and are accurately predicted by a comprehensive Langevin model. The diode-like behavior persists at the pore scale in geometries of higher dimension, where disorder impacts conductivity at the sample scale by extending the trapping regime over a broader range of flow speeds. This work has implications for our understanding of the survival strategies of magnetotactic bacteria in sediments and for developing their use in drug delivery applications in vascular networks.


Microfabrication
Corrugated microfluidic channels (height, h = 110 µm) are comprised of semi-circular walls with alternating concavity (Fig. 1). The microchannels were fabricated from polydimethylsiloxane (PDMS) using standard soft lithography methods 1 and plasma bonded to standard, 1 mm thick glass slides. 2 Walls forming the constriction at the channel throat have a fixed radius of curvature, R t = 35 µm, and a fixed wall separation width, w t = 70 µm. The radius of curvature of the expansion region, R p , which forms the pores, was varied, R p = [R t , 2R t , 3R t ], resulting in a pore width w p = w t + 2R t + 2R p = [3w t , 4w t , 5w t ] = [210 µm, 280 µm, 350 µm]. For 3D porous media experiments (Fig. 4), a suspension of polydisperse, transparent agarose hydrogel beads (radius, 10-25 µm; fine grade, Agarose Beads Technologies) was packed into a microfluidic channel (height, 75 µm; width, 1000 µm). The channel exit had a 45 • taper, which enabled the suspension to be jammed in the channel. After the packed bed was formed, the flow was switched to a suspension of magnetotactic bacteria under conditions described in the main text. Constant ambient flow maintained the arrested state of the bed for the duration of the experiment.

Microfluidics and pressure-driven flow calibration
To access the range of flow speeds relevant to the swimming magnetotactic bacteria within the microchannel throat, accurate control over the flow rate was required. A precision pressure controller (Elveflow OB-1) was used to drive the flow. The dynamic range of the applied pressure, P, was tuned by the addition of a serpentine resistor channel (height, 20 µm; width, 30 µm; total length, 10 cm) atop the corrugated test channel (Fig. S1). The resistor channels were plasma bonded to each test channel using standard methods. 2 We found that the resistor channels had a negligible effect on the bacterial image quality in the test channel.
For each experiment with magnetotactic bacteria, the physical control parameter for the cell dynamics is the maximal flow speed, U max , in the throat of the corrugated test channel. Thus, the maximal flow speed was calibrated as a function of applied pressure, P, for each test channel. Fluorescent tracer particles (0.5 µm diameter) in cell media were flowed through each corrugated channel at five different applied pressures. The flow field was reconstructed from video images of the tracer particles using custom particle tracking velocimetry (PTV) routines (MATLAB), and the maximal flow speed was measured in each throat (Fig. S3a-c, blue circles). The resulting maximal flow speeds were fitted as a function of pressure with the expected linear relationship, U max = a max P + b max , where the offset is due to hydrostatic pressure in the system (Fig. S3d-f, blue circles). The fitting coefficients are provided in Table S1. As a check, the same procedure was performed ( Fig. S3d-f, red squares) for the minimal flow speed, U min , along the channel centerline occurring in the pore (Fig. S3a-c, red squares). All experiments are performed in the low Reynolds number regime, ensuring linearity of the flow and pressure. In particular, the maximal Reynolds number of the bacterial experiments based on the channel height was Re = ρU max h/η ≤ 0.05, where ρ and η are the density and viscosity of water at room temperature, respectively. The ratio U max /U min was approximately constant to within ±2% for each geometry, consistent with the low Reynolds number flow conditions. For each experiment with cells, the corresponding U max was computed from this calibration. Finite element simulations (COMSOL) corroborate our PTV measurements of the flow topology ( Fig. S3g-i; see also below).

Langevin simulations in 1D porous media
Here, we illustrate the non-dimensionalization of the equations of motion for bacterial transport and discuss the key non-dimensional parameters that emerge, which are utilized throughout the main text. We then explain how the physical parameters characterizing magnetotactic motility are extracted from experimental data. We stress that there are no fitting parameters in the Langevin simulations, as all physical parameters are measured directly from experiments. Finally, we discuss the corresponding numerical simulation parameters.

Equations of motion for magnetotactic motility in flow
The Langevin equations of motion for magnetotactic bacteria in flow -originally introduced in the main text -are repeated here for convenience. This Langevin framework is well established for bacterial transport [3][4][5] and was recently adapted to magnetotactic bacteria. 6 The translational and rotational equations of motion in two dimensions arė Briefly, the orientation of the cells, θ , is measured relative to the upstream direction, x, and y is the cross-stream direction. The flow is described by the velocity field, u = (u, v), and the vorticity, ω = ∇ × u. Bacterial motility is parameterized by the mean cell swimming speed, V s , and rotational diffusivity, D r , where the latter characterizes fluctuations in cell orientation due to both Brownian motion and noise in flagellar motility. The cells have a rotational hydrodynamic mobility, µ r , and a magnetic moment, M, and they are exposed to a magnetic field, B = Be B , oriented in the positive x-direction.
The equations of motion are non-dimensionalized in terms of the cell swimming speed, V s , and the throat halfwidth of the channel, R 1D = w t /2, as follows: Here, the new dimensionless space and time variables are X = x/R 1D , Y = y/R 1D , and τ = tV s /R 1D . The linearity of the low Reynolds number flows considered here ensures that the velocity field topology does not change with the magnitude of the flow speed. The velocity field components and vorticity are thus scaled by the maximum value of the flow speed:ũ The key non-dimensional numbers that emerge from the system are where Φ is a non-dimensional flow speed, ∆ is a non-dimensional magnetic torque, andD r is a non-dimensional rotational diffusivity, or equivalently, an inverse rotational Péclet number.ξ is a normalized white Gaussian noise.

Boundary conditions
Steric boundary conditions prevent cells from entering the solid surfaces in our simulations. We implement a Weeks-Chandler-Andersen repulsive potential: where d is the distance to the closest solid surface, a is the bacteria radius. The hardness ε is set for its nondimensional valueε = ε/R 1D 6πηaV s = 2 × 10 −5 . The radius of the bacteria is set to be a = 0.5 µm, amounting to a non-dimensional valueã = a/R 1D = 0.014. This potential is zero away from the boundary and diverges fast as the cells gets closer to the surface. Cells are not oriented as the come close to the wall.

Experimental measurement of magnetotactic motility parameters
To ensure that the Langevin simulations are as physically representative of the experimental system as possible, all motility parameters of the swimming magnetotactic bacterial suspensions are measured directly, including D r , V s , and µ r M.

Rotational diffusivity
Micro-swimmers, including magnetotactic bacteria, experience a variety of noise in their swimming orientation, such as stochastic flagellar-actuated turning (e.g. run-and-tumble motility) and Brownian rotational diffusion. The resulting motility is often modeled as a continuous time random walk, where the noise is parameterized by an effective rotational diffusion, D r . The rotational diffusion coefficient is measured by tracking a suspension of free-swimming magnetotactic bacteria (439 trajectories) in the absence of an applied flow and magnetic field. The autocorrelation function of bacterial orientation, C ee (t) = e · e(t) , captures the persistence of the cell swimming direction, e, which decays with time (Fig. S2a). The decay or persistence time, τ p , for these cells was measured from an exponential fit to C ee (t) and found to be τ p = 4 s. The persistence time is directly related to the rotational diffusion coefficient as D r = 1/(2τ p ) = 0.125 s −1 .

Swimming speed
The mean swimming speed of the magnetotactic bacteria is also obtained by tracking a suspension of cells (> 1500 trajectories over 1 s long) in the absence of external flow, but with an applied magnetic field, comparable to the field strength used in the main text (B = 0.875 mT). The external magnetic field ensures that the main component of the swimming speed is co-planar with the imaging plane, and thus, there is no bias introduced by measuring a 2D projection of the 3D cell swimming velocity. The swimming speed distribution of the cells (Fig. S2b) had a mean value of V s = 135 µm/s and a standard deviation of 25 µm/s. The magnetotactic bacteria are robust swimmers, and we verified that their mean swimming speed decreases by only 20% over the course of one hour. Such variations in swimming speed have minimal impact on the relatively brief 30 minute experiments presented in the main text.

Magnetic alignment
Several techniques were previously established to estimate the magnetic moment, M, of magnetotactic bacteria. 7 These methods typically involve the measurement of the orientation of motile or non-motile cells as an imposed magnetic field is varied. Because the magnetic momentum is not perfectly aligned, trajectories will take an helical shape, 8 but their direction, when averaged over a few pitches, will follow the magnetic field. However, the coccoid MC-1 bacteria used in this study are spherical and small in diameter (d ≈ 1 µm), preventing direct optical measurement of the cell orientation. Furthermore, we also require the hydrodynamic rotational mobility, µ r , of the cells to solve the Langevin equations of motion. The mobility may be estimated analytically by assuming a spherical cell, but µ r ∼ d −3 potentially leading to large errors, and this approach ignores the hydrodynamic contribution of the flagella. Recognizing that the rotational mobility and magnetic moment only appear together in the equations of motion, their product, µ r M, is easily measured from controlled cell tracking experiments. We measure the ensemble average alignment of the swimming direction, e, with the magnetic field direction, e B , as e · e B for a range of magnetic field strengths, B. The result is predicted to obey the Langevin paramagnetic law: 7 where B 0 = D r /µ r M. The magnetic alignment of the cells is hence measured by tracking an ensemble of bacteria under a range of magnetic field strengths and computing the mean relative alignment, e · e B , which indeed follows the Langevin paramagnetic law (Fig. S2c). The result is fitted by the Langevin function to extract the characteristic magnetic field B 0 = 45 µT, which is relatively close to Earth's magnetic field of 56 µT. Because the rotational diffusivity of the cells is measured independently (see above), we can obtain an estimate for the product of the magnetic moment and rotational hydrodynamic mobility as µ r M = D r /B 0 = 2800 T −1 s −1 .

Simulation parameters
To demonstrate this first observation of vortical localization of upstream swimming cells in porous media flows, we vary both the geometry and the flow strength, while holding the magnetic field fixed (B = 1.75 mT). The relevant non-dimensional magnetic alignment parameter is thus ∆ = R 1D MBµ r /V s = 1.25, and the non-dimensional rotational diffusion coefficient isD r = D r R 1D /V s = 0.032 for all simulations. The non-dimensional flow strength is varied in the range Φ = U max /V s ∈ [0, 10], and the geometries tested had pore to throat aspect ratios corresponding to experiments, including w p /w t = [3,4,5], as well as straight channels. Experimental particle tracking velocimentry (PTV) measurements of the channel flow fields (Fig. S3a-c) are used primarily to calibrate the maximal flow strength, U max , as a function of applied pressure, P, (Fig. S3d-f) and to validate finite element simulations (COMSOL) of the flow fields ( Fig. S3g-i). For Langevin simulations of the bacterial dynamics, we rely on the finite element simulations of the flow fields in each geometry. These two-dimensional, finite element simulations are a good approximation to the experimental flow fields, especially at the high-aspect ratio throat, and provide higher resolution of u(x, y), v(x, y), and ω(x, y). The finite element simulations are computed in the 2D channel domain using Stokes equation, 0 = −∇p + η∇ 2 u, and continuity, ∇ · u = 0, with a no slip boundary condition on the corrugated walls for a single, arbitrary applied pressure gradient for each geometry. The finite element simulation domain consisted of 11 spatial corrugation periods with prescribed constant pressure conditions at the inlet and outlet, and it was meshed with over 70,000 vertices in each geometry. The velocity field from the most central corrugated pore is used as the periodic flow field for the Langevin simulations, and the vorticity, ω, is computed directly from this flow field. Linearity enables the simulated flow field and vorticity field topologies to be scaled with the maximal centerline flow speed, U max , for each Φ tested.
For the Langevin simulations, the non-dimensional time step was fixed at δ τ = 0.005. Time integration was performed by a fourth-order stochastic Runge-Kutta scheme, and the flow and vorticity fields were interpolated at each time step by a bilinear interpolation. Periodic boundary conditions were implemented in the stream-wise direction. For each condition, we simulated 600 cell trajectories for a total dimensionless simulation time, τ ≥ 20,000.

Influence of the magnetic field orientation on the bacterial conductivity
In our experiments, we used a magnetic field direction oriented anti-parallel to the flow. While this idealized system yields striking observations of dynamical cell trapping, in general, the flow direction may have an arbitrary orientation relative to the field, for example in environmental conditions or biomedical applications. We thus exploit Langevin simulations to determine the effect of the magnetic field orientation relative to the flow on the trapping phenomena demonstrated (Fig. S7) in the main text. In the corrugated channels, without flow, cells are directed by the magnetic field and guided by mechanical interactions with walls ( Fig. S7f-j). When cells are directed into a concave wall of the corrugated channel, their random walk motility relative to the field strength is generally insufficient to enable escape from the wall. However, this idealized geometry was specifically designed to observe the vortical trapping effect illustrated in the main text ( Fig. 1), which can be clearly observed for moderate flows (Fig. S7k-o). The bias in cell swimming direction creates an asymmetry in the vortical structure, which only appears on one side of the throat. By driving the cells on the side of the corrugated channels where the fluid velocity is also weaker, we observe that the plateau of the trapping regime is significantly extended (Fig. S7a-e), and larger flow rates are needed to wash away the cells and enter the downstream regime ( Fig. S7p-t). In conclusion, the trapping regime is increased by the introduction of a relative angle between the magnetic field and flow direction, it is because of a cooperative effect between (i) the vortical trapping, which still occurs in these tilted configurations, (ii) steric trapping where cells are driven against the concave boundary of the channel (reducing the upstream regime), and (iii) their localization in a low flow speed region of the channel (reducing the downstream regime; Fig. S7a-e).

Langevin simulations in 2D porous media and the effect of disorder 2.3.1 Geometry, flow, and Langevin simulations
In this section, we use the numerical Langevin model detailed above to study in detail the the effects of disordered porous media in higher spatial dimension on cell transport for swimming magnetotactic bacteria directed upstream against a fluid flow.
The geometries considered here are based on an array of circular obstacles 5 (Fig. S8a). We consider both an ordered and a disordered arrangement, where the ordered geometry comprises a 10 × 10 square lattice of cylinders with porosity 0.5. The cylinders have radius, R 2D = 105 µm, which is three times the half-width of the throat for the 1D corrugated system, R 1D , with a center to center distance a 2D = 8R 1D = 280 µm. The wall-to-wall distance between two neighboring obstacles is hence equal to the throat width of the 1D corrugated medium, 2R 1D = 70 µm.
The disordered medium comprises the same number of obstacles and same porosity as the ordered medium, but obstacles are randomly distributed (non-overlapping; Fig. S8b).
For the purposes of simulations, the disordered and ordered media are made globally periodic by tiling the obstacle geometry in a 3 × 3 array. Flow fields are computed numerically for each medium using Stokes flow assumptions in COMSOL as detailed above. The 3×3 array was placed in the middle of a rectangular channel with a length of 161.2 mm and width of 81.2 mm leaving an open space with a length of 38.95 mm both at the inlet and outlet of the channel. The maximum Reynolds number in the 2D porous geometries was Re ≤ 0.1, which suggests that the entrance length was L ent = 0.05 LRe ≤ 0.35 µm for the highest flow speed examined, where L is taken as the wall-to-wall distance between two neighboring obstacles in the ordered lattice. The entrance length was thus sufficiently smaller than the length of the open space at the inlet as well as the length of a single tiling of the pillar arrangement (i.e. each one of the 3×3 array) ensuring a fully developed flow profile at the central portion of the 3×3 array. The simulations used a combination of triangular and quadrilateral mesh elements (2,605,814 elements for ordered and 4,111,650 elements for disordered geometry), which was greatly refined near the pillar surfaces (≈ 2 µm) to resolve the high velocity gradients. The upstream boundary condition was a constant unidirectional flow imposed at the channel inlet, and the downstream boundary condition at the outlet was a reference value for the pressure set to zero. No-slip boundary conditions were implemented on the pillar surfaces and slip boundary conditions on the sidewalls of the channel.
The central tile of the 3×3 flow field was used for Langevin simulations with periodic boundary conditions for both the fluid flow and simulated cell dynamics. For Langevin simulations, 2500 cells are initially placed homogeneously in the porous media in the disordered geometry, (400 cells in the ordered geometry). Their motion is then simulated via the Langevin equations above using the same Runge-Kutta scheme, same wall/obstacle boundary conditions, same time steps, and the same (non-dimensional) total simulation time τ = tV S /R 1D = 256 for the disordered geometry, and τ = tV S /R 1D = 64 in the ordered geometry (See 9). The sample scale conductivity at long times is defined from the ensemble speed average V (t f ) at the final time t f as V (t f ) /V S .

Population dynamics reveal that heterogeneous flow in disordered media enhances trapping at the sample scale
Due to the geometrical uniformity of the flow field in the ordered media, all pores/throats uniformly exist in one of the three identified transport regimes (Fig. 4, Fig. S9a-b) corresponding to the 1D case (Fig. 2). However, the heterogeneity of the flow field in the disordered geometry enables the coexistence of the upstream, downstream, and trapping regimes (Fig. 4, Fig. S9c-d). To elucidate the role of disorder, we take a cell population dynamics approach. We associate each simulated cell with an instantaneous transport state -upstream, trapped, or downstream -depending on its velocity over 3R 2D /V s . Cells with a net displacement ∆X < −R 2D /2 during this time are classified as in the downstream regime, while cells with a displacement ∆X > R 2D /2 are classified as in the upstream regime. In between, cells are classified as being in the trapped regime, −R 2D /2 ≥ ∆X ≥ R 2D /2. The evolution of these cell populations from an initially spatially uniform distribution is then studied during the simulations (Fig. S9, top row). We represented the initial and final spatial distributions of the cells and their states for three different dominant sample scale conductivity regimes for both geometries, i.e., V /V S (t ∞ ) = 0.5, 0, and -0.5 (Fig. S9, rows 2-4).
In their initially uniformly seeded positions, cells coexist in three corresponding populations over a wide range of flow speeds in the disordered geometry (Fig. S9c, top row), when compared to the ordered geometry (Fig. S9a, top  row). In the ordered geometry, individual cells quickly settle into the regime associated with the sample scale due to the uniformity of the flow topology across pores, for example the trapping regime shown in Fig. S9a-b (row 3). In contrast, the state of the cell populations in the disordered geometry evolve more slowly: As cells circulate and sample the heterogeneous flow field (Fig. 4), they have a propensity to become trapped in specific pores, for example Fig. S9cd (row 3). The proportion of trapped cells increases over time for all flow regimes investigated (Fig. S9c-d, rows 2-4). This feature emphasizes the fact that the residency time within the trapping regions becomes very long compared to any other time scales involved in cell transport, including upstream or downstream excursions. Cells circulate fast through these pores in the upstream and downstream states and conversely, reside for an extended time in trapping pores. Even though the trapping is a pore scale effect, the spatial correlations of the flow field (larger than the pore scale) translate into patches of pores in the same regime (Fig. S9c, third row, blue patch). Thus, while the nonlinear conductivity is a pore scale effect, the correlated heterogeneity of porous media flows results in spatial correlation of the bacterial conductivity. Tables   Supplementary Table 1: Fitting coefficients for the maximal, U max , and minimal, U min , flow speed versus pressure drop curves (Fig. S3d-f) Figure 1: Schematic of the corrugated microfluidic device and resistor channel. A serpentine resistor channel (height, 20 µm; width, 30 µm; total length, 10 cm) regulates the dynamic range of the applied pressure, and thus, controls the flow speed within the corrugated test channel. The resistor and corrugated test channels are both fabricated from PDMS. The outlet of the resistor channel is aligned with the inlet of the corrugated channel, and the two PDMS structures plasma bonded together. The resistor and corrugated test channels are also bonded to a standard glass microscope slide (thickness, 1 mm; not to scale).  , For each of the three corrugated microchannels, the flow strength was calibrated with the applied pressure for both the maximum and minimum centerline flow speeds. A linear fit to five experimental measurements, U min,max = a min,max P + b min,max , provides the calibration coefficients, which are listed in Table S1. g-i, Finite element simulations (COMSOL) of the flow field for each of the three corrugated channels were performed in the low Reynolds number Stokes flow limit. The simulated flow fields corroborate PTV measurements and are utilized in subsequent Langevin simulations. Scale bars, 35 µm.   Figure 7: Influence of magnetic field angle with the flow direction on pore scale transport. a-e, Langevin simulations reveal a systematic broadening of the trapping regime for the mean cell speed as the magnetic field angle increases relative to the flow direction for angles 0 • , 15 • , 30 • , 45 • , 75 • . Green, red, and blue curves correspond to large, medium, and small corrugations (see main text). Vertical black lines mark the upstream (Φ 1 ), trapping (Φ 2 ), and downstream (Φ 3 ) regimes observed for 0 • in the medium corrugation (see also Fig. 2). f-t, For each three regimes defined in the main text, trajectories are represented for the medium corrugation corresponding to a-e.  Figure 9: Disorder promotes the coexistence of upstream, trapping, and downstream transport regimes a-d, (row 1) Populations of cells in the upstream (green), trapping (red), and downstream (blue) regimes. Initial (a, t 0 ) and final (b, t ∞ ) conditions for the ordered geometry and disordered geometry (c-d) are shown, respectively. Upper colored ticks indicate the average cell flow V /V s rates value of 0.5, 0 and -0.5, colors corresponding to main text (Fig. 4-h). a-d, (row 2) Evolution of the cell population and cell locations in the upstream regime for the ordered geometry V (t ∞ )/V s = 0.5 (a-b) and for the disordered geometry (c-d). a-d, (row 3) Same as row 2 for the trapping regime V /V s = 0 as defined by the ordered geometry. a-d, (row 4) Same as rows 2-3 for the downstream regime V /V s = −0.5 as defined by the ordered geometry.