Microbial transport and dispersion in heterogeneous flows created by pillar arrays

We modeled a microbe as a prolate ellipsoid which has an effective aspect ratio $q$, a constant swimming speed $V$, and an effective rotational diffusion coefficient $D_{r}$. In the 2D space, the flow field at position $(x,y)$ is denoted as $\bm{U_{f}}(x,y)=[u(x,y),\,v(x,y)]$, where $u$ is the flow speed along $x$ direction and $v$ is the flow speed along $y$ direction. The microbe performs a run-and-tumble motion, determined by the flow speed $\bm{U_{f}}$, the swimming speed $V$, gradients of the flow field (i.e. $\dfrac{\partial u}{\partial x}$, $\dfrac{\partial u}{\partial y}$, $\dfrac{\partial v}{\partial x}$, $\dfrac{\partial v}{\partial y}$) and the rotational diffusivity $D_{r}$. The equations of microbial motion in the 2D flow are expressed as [17]

Where $(x,y)$ is the microbial position. $\theta$ is the microbial moving direction, which is also the direction of microbial major body axis. $u$ is the background flow speed along $x$ direction. $v$ is the background flow speed along y direction. $V$ is the microbial swimming speed. $q$ is the aspect ratio. And $\xi_{r}$ is the rotational noise represented as a Gaussian-distributed angular velocity with mean zero and variance $2D_{r}/\Delta t$, $D_{r}$ is the rotational diffusion coefficient.

The equations of motion were solved using a fourth-order Runge-Kutta scheme implemented in Python. Trajectories of $10^{7}$ microbes are simulated. Based on the experimental measurements [17], for P. aeruginosa, the microbial parameters are set as $q=9.4$, $V=45\,\mu m/s$ and $D_{r}=1.4\,rad^{2}/s$. The time step $\Delta t$ is chosen as 0.001 $s$. The space discretization $\Delta x$ is chosen as 0.05 $\mu m$.

The background flow field $\bm{U_{f}}=[u,v]$ influenced by the size and arrangement of pillars is computed by numerically solving the 2D stokes equation, $\mu\nabla^{2}\bm{U_{f}}-\nabla p=0$, in COMSOL Multiphysics, where $\bm{U_{f}}$ is the fluid velocity field, $p$ is the pressure and $\mu$ is the dynamic viscosity of water. Specifically, fluid flows through a rectangular chamber (500 $\mu m\times 1000\,\mu m$) whose size is large enough such that the chamber boundaries have negligible effect on the microbial movements in the simulation time. Inside the chamber, there is an array of cylindrical pillars (5 columns, 10 rows). The pillar size is characterized by its radius $R$. Adjacent rows and columns of pillars are separated by a distance $W$ ($W=100\,\mu m$) in both $x$ and $y$ directions, as shown in Figure 1a. In addition, in an asymmetric array of pillars, adjacent rows shift a distance $dH$ in the $y$ direction. Fluid flows with a constant speed from the left side of the chamber into the device. The right side of the chamber is set as a stress-free condition. And the two other sides of the chamber are set as periodic boundary conditions. And the pillar surfaces are set as the no-slip boundaries. The microbes flow into the pillar array from the center of the left side of the chamber where 200 $\mu m\leq y\leq 300\,\mu m$.

The no-slip pillar surface creates local velocity gradients of the flow and influences microbial trajectories. To quantify the flow properties, distributions of shear rate and Lagrangian fluid stretching are calculated. The shear rate $\dot{\gamma}$, representing the velocity gradients, is calculated as the positive eigenvalues of the strain rate tensor $\bm{E}=\dfrac{1}{2}(\nabla\bm{U_{f}}+\nabla\bm{U_{f}}^{T})$.

The Lagrangian fluid stretching, which quantifies the deformation of a small circular fluid element into an ellipse due to the advection of ambient flow field, is calculated as the square root of the maximum eigenvalue of the left Gauchy-Green deformation tensors. Specifically, suppose a fluid particle which locates at position $\bm{x_{0}}$ at time $t=0$ is deformed by the flow field $\bm{U_{f}}$ and advected to a new position $\bm{\Phi}_{\lambda}(\bm{x_{0}})$ at time $t=\lambda$, where $\bm{\Phi}_{\lambda}=\bm{x_{0}}+\int_{0}^{\lambda}\bm{U_{f}}dt$ is the flow map. The Lagrangian fluid stretching is thus calculated as the eigenvalues of $\bm{C}_{\lambda}(\bm{x_{0}})=(\nabla\Phi_{\lambda}(\bm{x_{0}}))(\nabla\Phi_{\lambda}(\bm{x_{0}}))^{T}$. In a practical manner, we followed the approach of Parsa et al. [31] to compute the Lagrangian stretching field. First, a regular grid of virtual particle position $\bm{x}$ ($500\times 500$) within a unit cell of the pillar array (which contains one pillar) is chosen. For each virtual particle, we define four auxiliary points around it. The auxiliary points are located a small distance (0.01 $\mu m$) to the north, south, east and west from the virtual particle. By numerical integration of the flow field $\bm{U_{f}}$ backward in time (i.e. $-\lambda$), the trajectories of the virtual particle and the auxiliary points are obtained to get the flow map $\bm{\Phi_{\lambda}}$. During the numerical integration, every 0.01 $s$, the trajectories are rescaled to ensure the distances between auxiliary points and the virtual particle are always 0.01 $\mu m$. The flow map gradients, $\bm{\Phi}_{\lambda}(\bm{x_{0}})$, are then computed via the central differences using the four auxiliary points surrounding $\bm{x_{0}}$.

To quantify the microscopic microbial dynamics, the rotational Péclet number, defined as $Pe_{r}=\dfrac{<|\omega|>}{2D_{r}}$, is calculated, where $<|\omega|>$ is the mean absolute vorticity ($\omega=\nabla\times\bm{U_{f}}$), and $D_{r}$ is the rotational diffusion coefficient of the microbes. After obtaining the microbial trajectories by solving equations (1)-(3), the mean square displacement (MSD), the effective transverse velocity (the velocity perpendicular to the flow), the effective transverse dispersion coefficient, the effective lateral velocity (the velocity parallel to the flow) and the effective lateral dispersion coefficient are determined. The MSD is calculated as

where $N$ is the number of microbes ($N=10^{7}$ in this study), $(x_{i}(t),y_{i}(t))$ is the position of the microbe $i$ at time $t$.

The maximum likelihood estimator (MLE) is used to infer the velocities and dispersion coefficients from microbial trajectories [32, 33]. Specifically, suppose the microbial distribution along $x$ direction follows $p(x_{i};\,v_{x},\,D_{x})=\dfrac{1}{\sqrt{4\pi D_{x}\delta}}e^{-\frac{(x_{i}-v_{x}\delta)^{2}}{4D_{x}\delta}}$, where $x_{i}$ is the displacement of microbe $i$ along $x$ direction (i.e. $x_{i}=x_{i}(\delta)-x_{i}(0)$), $\delta$ is the time the microbe takes to achieve displacement $x_{i}$, $v_{x}$ is the lateral velocity, and $D_{x}$ is the lateral dispersion coefficient. The likelihood for obtaining $N$ microbial displacements $\{x_{i}\}$ is expressed as

By solving $\arg\max f(x_{i};\,v_{x},\,D_{x})$, we get the estimated lateral velocity $\hat{V}_{x}$ and the estimated lateral dispersion coefficient $\hat{D}_{x}$ as

Similarly, based on the microbial displacements in the $y$ direction ${y_{i}}$, the transverse velocity $\hat{V}_{y}$ and the transverse dispersion coefficient $\hat{D}_{y}$ can be estimated using the MLE as

We investigated the effect of heterogeneous flow created by the pillar array on swimming microbes by tracking their trajectories and distributions. In a rectangular chamber consisting of a periodic array of circular pillars, fluid flows from the left side of the chamber to the right, and microbes enter the chamber from the center of the left inlet and swim with the flow (Figure 1a). The pillar radius $R$ varies from 10 $\mu m$ to 40 $\mu m$. The spacing between two adjacent pillars $W$, which is also the width and height of one periodic unit, is held constant ($W=100\,\mu m$). The mean fluid flow speed varies from 0 to $400\,\mu m/s$, corresponding to a mean shear rate $\dot{\gamma}_{m}$ approximately ranging from 0 to 10 $s^{-1}$. Figure 1b shows the flow streamlines in one unit. Due to the symmetry, the position a streamline enters the unit, denoted as $y_{i}$, is the same as the position it leaves the unit, denoted as $y_{i+1}$. Therefore, small-size passive particles, including the non-spherical particles, flow through the unit following the streamlines. And the particle distribution at the right outlet of the unit is the same as the distribution at the left inlet when diffusion is negligible. However, since the ellipsoidal microbes have active swimming ability, their trajectories are deviated from the flow streamlines ($y_{i+1}\neq y_{i}$). Specifically, the fluid shear $\dot{\gamma}$, as shown in Figure 1c and f, align the major axis of the microbial elongated body perpendicular to the pillar surface. In the vicinity of the pillar, the swimming speed of the microbe is larger than the local flow velocity, so the microbes will be able to swim upstream and reach the leeward side of the pillar. Moreover, to consider the nonlocal historical effects stemming from advection, we further calculated the Lagrangia fluid stretching field over a time interval $\lambda=3.5\,s$ (Figure 1d and g). Since the microbial density is strongly correlated with the regions of high stretching [26], Figure 1d and g indicate that microbes that locate close to the pillar tend to be attracted to swim around the pillar. When the pillar radius increases, the influential distance within which the microbes might be captured by the pillar becomes larger. Figure 1e and h show that a large number of microbes accumulate behind the pillar due to the bending of microbial trajectories by the local fluid shear. And the time for most microbes arriving at the leeward side of the pillar is influenced by its radius $R$. Generally, when the pillar radius increases, the population of microbes move slower around the pillar, and the time microbes take to accumulate at the leeward pillar side is larger (Figure S1). However, by tracking the position where a microbe passes line AA’ (line AA’ locates between the first and second columns of pillars), we found that the pillar radius has only slight influence on the microbial distribution (Figure 1e and h, Figure S1). When $R=0.1W$ and $R=0.25W$, microbial distributions have a large peak at the center and heavy tails at the two sides (Figure 1e and h). In addition, microbes that locate far from the pillar pass line AA’ earlier than microbes near the pillar surface (Figure S1). In a macroscopic point of view, the reduction of the flowing speed of microbial population by the pillar geometries is reminiscent to the dynamics of highly viscous fluid flowing through the pillar array. Nevertheless, the fluid viscosity arises from the internal frictional force between adjacent layers of fluid. But the slower changes of microbial densification patterns are caused by the interactions between individual microbes and the pillar surface.

We further quantified the effects of pillar radius and fluid mean shear rate on the microbial motility. Specifically, the MSD of a particle can be approximated by the power law as $MSD=D_{\alpha}t^{\alpha}$, where $D_{\alpha}$ is the equivalent diffusion coefficient, and $\alpha$ is the exponent[34]. Our results show that $\alpha$ is close to 2 regardless of the radius $R$ and the mean fluid shear rate $\dot{\gamma}_{m}$. $D_{\alpha}$ increases over $\dot{\gamma}_{m}$ and decreases over $R$ (Figure 2a and b). Particularly, a maximum likelihood estimator (MLE) is applied to analyze the advection speeds and dispersion coefficients of microbes. We found that the estimated lateral speed $\hat{V}_{x}$ and the lateral dispersion coefficient $\hat{D}_{x}$ increases when fluid flows faster in the pillar array (i.e., when $\dot{\gamma}_{m}$ increases), but the transverse dispersion coefficient $\hat{D}_{y}$ decreases over $\dot{\gamma}_{m}$. This result agrees well with previous findings that hydrodynamic gradients hinder transverse bacterial dispersion in a microfluidic crystal lattice [26]. Moreover, we found that the pillar radius $R$ has a larger influence on the lateral microbial dynamics ($\hat{V}_{x}$ and $\hat{D}_{x}$) when $\dot{\gamma}_{m}$ is large, but influences transverse dynamics $\hat{D}_{y}$ significantly when $\dot{\gamma}_{m}$ is small. Large-size pillars hinder both lateral and transverse microbial movements, especially $\hat{V}_{x}$ and $\hat{D}_{y}$. When the pillar has an intermediate radius and $\dot{\gamma}_{m}$ is large, the lateral dispersion coefficient $\hat{D}_{x}$ is the largest (Figure 2c, d and e).

To reveal the mechanisms that cause the transient microbial accumulation and the slow changing of microbial density surrounding the pillar, the microscopic dynamics of microbes is analyzed. The trajectories of microbes in the pillar arrays are first visualized. Figure 3a shows that when the pillar radius $R$ is small, many microbes maintain relatively straight trajectories and flow through the array quickly. Only some trajectories that are close to the pillar are redirected to move around the pillar. However, when the pillar radius $R$ increases, meandering trajectories will dominate. Specifically, three types of movements, ‘swinging’, ‘zigzag’ and ‘adhesive’ movements are observed. When a microbe locates between two rows of pillars, it tends to swing up and down while keeping moving forward. This microbial motion is thus referred to as the ‘swinging’ motion (Figure 3b (1)). When a microbe reaches the vicinity of the pillar, the local velocity gradient applies hydrodynamic torque on the microbe, leading to the frequent loops near the pillar surface (Figure 3a). Though the microbe keeps rotating, it is hardly to swim far away from the top or the bottom sides of the pillar where the shear rate is the largest, as shown in Figure 1f. If the microbe continues swimming from the bottom of the pillar to its leeward side and shifts to the top side of the next pillar, this type of movement is referred to as a ‘zigzag’ mode (Figure 3b (2)). If the microbe moves along the top pillar surface and shifts to the top side of the next pillar, it is then referred to as an ‘adhesive’ mode (Figure 3b (3)). To understand how the pillar radius $R$ and the mean shear rate $\dot{\gamma}_{m}$ influence the three types of movements, the microbial moving direction $\phi$ in a circular coordinate, where the coordinate center is the pillar center, is calculated. When $\dot{\gamma}_{m}=6\,s^{-1}$, there is a thin dark layer surrounding the pillar surface, where $\phi$ is close to $\pi/2$ (Figure 3c). In this layer, microbes move parallel to the pillar surface, resulting in a ‘zigzag’ or an ‘adhesive’ movement. When $R=0.3W$ and $\dot{\gamma}_{m}=6\,s^{-1}$, there is a white area next to the thin dark layer (Figure 3c). In this white area, microbes are moving away from the pillar, which is likely to lead to a ‘swinging’ motion. The white area is hard to be recognized when $R=0.1W$ (Figure 3c). Thus, the ‘swinging’ motion is less frequent when pillar size is reduced. When $\dot{\gamma}_{m}=1\,s^{-1}$, the dark thin layer and the white area disappear (Figure 3c), suggesting the microbial trajectories become more randomized. Moreover, the cell density along the red line in Figure 3d in the $y-\theta$ phase space is calculated. When $R=0.3W$ and $\dot{\gamma}_{m}=6\,s^{-1}$ (Figure 3e, left), around half of the microbes move toward the pillar leeward side, which is represented by the colorful areas in the left bottom and right upper corners. When the radius R decreases to $0.1W$, the portion of microbes that flow horizontally becomes larger (Figure 3e, right). And when $\dot{\gamma}_{m}=1\,s^{-1}$, the two probabilities of a microbe moving toward and away from the pillar get closer (Figure 3e, central). The meandering trajectories and the frequent loops increase the distance a microbe travels, such that the speed of microbial population moving forward is reduced. Figure 3f shows that a larger radius $R$ and a small shear rate $\dot{\gamma}_{m}$ lead to long moving distances of microbes to reach the areas behind the pillar. Besides, the slower speed of the microbe near the pillar surface (Figure 3a) also affects the speed of the microbial population. Therefore, in the array of large-size pillar, microbial transient accumulation and the slow change of microbial density is observed (Figure S1).

When the geometric symmetry is broken, few previous experiments [21, 26, 35, 36] has shown that the microbial transport and pattern are influenced by the asymmetric structures. There are two general ways to break the symmetry of a pillar array. One is to rotate the incidence angle of the flow relative to the pillar array, which is also referred to as the rotated square lattice. In this scenario, Dehkharghani et al. [26] found that the microbial densification patterns and the dispersion coefficients are affected by the incidence angle. Another is to shift adjacent columns of pillars by a small distance $dH$ and create a stretched array, which is referred to as a stretched DLD array (or a row-shifted parallelogram array). The stretched DLD array has been used to separate rod-like bacteria [22, 28, 29]. However, since the results in the rotated square lattice layout do not extend to the stretched DLD array, the dynamics and regulatory mechanisms of microbes in the stretched DLD array are still missing.

As shown in Figure 4a, the stretched DLD pillar array contains rows of pillars with radius $R$. Adjacent rows and adjacent columns are separated by a distance $W$. The right neighboring pillar is shifted a distance $dH$ in the $y$ direction. The shift between one column and its $M$ nearest columns is $MdH$, which is chosen to coincide with $W$. Thus, the array is periodic and invariant to a translation $MdH$ in the $y$ direction. And the stretched DLD array has a built-in angle $\theta_{DLD}=\arctan{1/M}$. Fluid flows in the horizontal direction into the pillar array, resulting in a large flow speed between two rows of pillars and a small speed behind each individual pillar (Figure 4b). Two adjacent ellipsoidal areas that have the largest flow velocity $u$ is also shifted $dH$ in the $y$ direction. Microbial trajectories bend around the pillar and are influenced by the built-in angle $\theta_{DLD}$ (Figure 4c). Compared with the symmetric pillar array, the distribution of the shear rate $\dot{\gamma}$ is rotated around the pillar center by $\theta_{DLD}$ (Figure 4d, g and j). The distribution of Lagrangia stretching becomes asymmetric about the $x$ axis (Figure 4e, h and k). When $R=0.3W$ and $\dot{\gamma}_{m}=6\,s^{-1}$, the large-stretching regions at the top side of the pillar extend away from the pillar surface, while at the bottom side they concentrate near the pillar surface (Figure 4e). The asymmetric fluid properties induce the macroscopic movements of microbes in the $y$ direction. Figure 4f shows that after passing four columns of pillars, the peak of the microbial density shifts toward the bottom. When $R$ decreases to $0.1W$, the Lagrangia stretching is still large surrounding the pillar surface. But at the top side of the pillar, two streamline-like regions that have large Lagrangia stretching appear. At the bottom side, the Lagrangia stretching remains small (Figure 4k). Consequently, microbes that have the highest density at the center start to separate into two groups. Around 45% of microbes are in the upper group and move faster than the rest 55% of microbes in the bottom group. When the mean shear rate $\dot{\gamma}_{m}$ decreases (i.e. $\dot{\gamma}_{m}=1\,s^{-1}$), the influence of the asymmetric pillar arrangement on the flow field is small (Figure 4h), and the microbial movements are more randomized (Figure 4i).

To get a quantitative understanding of the microbial motility in the stretched DLD pillar array, we then estimate the lateral and transverse microbial velocities and dispersion coefficients based on a maximum likelihood estimator (MLE). Results in Figure 5 show that the velocities and the dispersion coefficients are either linearly or parabolically proportional to the parameters that modulate the heterogeneous flow field, including the pillar radius $R$, the asymmetric shift $dH$ and the mean shear rate $\dot{\gamma}_{m}$. Specifically, when the pillar radius increases, the lateral velocity along the flow direction decreases in a parabolic manner (Figure 5a), while the transverse velocity increases linearly (Figure 5b). The dispersion coefficients have the opposite trend over the pillar radius. The lateral dispersion increases (Figure 5c), but the transverse dispersion is reduced over the radius (Figure 5d). When the asymmetric shift becomes large, the lateral velocity and the dispersion coefficient decreases parabolically (Figure 5e and g), while the transverse motions ($\hat{V}_{y}$ and $\hat{D}_{y}$) are enhanced (Figure 5f and h). As for the mean shear rate, when fluid flows faster ($\dot{\gamma}_{m}$ increases), the lateral and transverse velocities increase linearly (Figure 5i and j). The lateral dispersion coefficient also increases, but in a parabolic manner (Figure 5k). The transverse dispersion is reduced (Figure 5l). Overall, these results indicate that the structural asymmetry enhances the transverse microbial movements ($\hat{V}_{y}$ and $\hat{D}_{y}$). The flow speed (reflected by $\dot{\gamma}_{m}$) increases microbial advection ($\hat{V}_{x}$ and $\hat{V}_{y}$). And large pillars hinder the lateral advection and the transverse diffusion, but improve the transverse advection and the lateral diffusion.

Since the lateral and transverse movements of microbes are the macroscopic reflections of interactions among microbial swimming motion, fluid flow and structural boundaries, the microbial properties should also play significant roles to determine its advection and dispersion. Therefore, we further quantify the effects of three microbial properties, including the aspect ratio $q$, which describes how prolate the microbial body is, microbial swimming speed $V$, and the rotational diffusion coefficient of the microbe $D_{r}$. Results in Figure 6a and b show that prolate microbes tend to have larger advection speed in both $x$ and $y$ directions. The lateral diffusion coefficient decreases parabolically over the aspect ratio $q$ (Figure 6c). The transverse diffusion coefficient is less influenced when $q$ is larger than 4 (Figure 6d). The swimming speed $V$ has linear effects on the microbial motility. Surprisingly, velocities $\hat{V}_{x}$, $\hat{V}_{y}$ and the lateral diffusion coefficient $\hat{D}_{x}$ decreases over the swimming speed $V$ (Figure 6e, f and g). Only the transverse diffusion coefficient $\hat{D}_{y}$ increases over $V$ (Figure 6h). This result suggests that the fast moving of microbial population in the pillar array is mainly attributed to the flow instead of the swimming of microbe individuals. In addition, the lateral microbial movements ($\hat{V}_{x}$ and $\hat{D}_{x}$) increase when the rotation (or tumbling) of the microbe is enhanced (Figure 6i and k). The transverse velocity $\hat{V}_{y}$ is less influenced when the rotational diffusion coefficient $D_{r}$ is large (Figure 6j). And the transverse dispersion $\hat{D}_{y}$ decreases over $D_{r}$ (Figure 6l). It should also be noted that the influence of the aspect ratio $q$ is much smaller than the influence of the swimming speed $V$ and the rotational diffusion coefficient $D_{r}$ (Figure 6). In summary, prolate microbes that swim slowly but can tumble quickly are likely to have large advection velocities and also large lateral dispersion in the stretched DLD pillar arrays. In contrast, prolate microbes that swim fast and tumble less frequently usually have large transverse dispersion but small velocities and small lateral dispersion.