Publications

It is well-known that by placing judiciously chosen image point forces and doublets to the Stokeslet above a flat wall, the no-slip boundary condition can be conveniently imposed on the wall [Blake, J. R. Math. Proc. Camb. Philos. Soc. 70(2), 1971: 303.]. However, to further impose periodic boundary conditions on directions parallel to the wall usually involves tedious derivations because single or double periodicity in Stokes flow may require the periodic unit to have no net force, which is not satisfied by the well-known image system. In this work we present a force-neutral image system. This neutrality allows us to represent the Stokes image system in a universal formulation for non-periodic, singly periodic and doubly periodic geometries. This formulation enables the black-box style usage of fast kernel summation methods. We demonstrate the efficiency and accuracy of this new image method with the periodic kernel independent fast multipole method in both non-periodic and doubly periodic geometries. We then extend this new image system to other widely used Stokes fundamental solutions, including the Laplacian of the Stokeslet and the Rotne-Prager-Yamakawa tensor.

We show using theory and experiments that a small particle moving along an elastic membrane through a viscous fluid is repelled from the membrane due to hydro-elastic forces. The viscous stress field produces an elastic disturbance leading to particle-wave coupling. We derive an analytic expression for the particle trajectory in the lubrication limit, bypassing the construction of the detailed velocity and pressure fields. The normal force is quadratic in the parallel speed, and is a function of the tension and bending resistance of the membrane. Experimentally, we measure the normal displacement of spheres sedimenting along an elastic membrane and find quantitative agreement with the theoretical predictions with no fitting parameters. We experimentally demonstrate the effect to be strong enough for particle separation and sorting. We discuss the significance of these results for bio-membranes and propose our model for membrane elasticity measurements.

Cilia and flagella are essential building blocks for biological fluid transport and locomotion at the micrometre scale. They often beat in synchrony and may transition between different synchronization modes in the same cell type. Here, we investigate the behaviour of elastic microfilaments, protruding from a surface and driven at their base by a configuration-dependent torque. We consider full hydrodynamic interactions among and within filaments and no slip at the surface. Isolated filaments exhibit periodic deformations, with increasing waviness and frequency as the magnitude of the driving torque increases. Two nearby but independently driven filaments synchronize their beating in-phase or anti-phase. This synchrony arises autonomously via the interplay between hydrodynamic coupling and filament elasticity. Importantly, in-phase and anti-phase synchronization modes are bistable and coexist for a range of driving torques and separation distances. These findings are consistent with experimental observations of in-phase and anti-phase synchronization in pairs of cilia and flagella and could have important implications on understanding the biophysical mechanisms underlying transitions between multiple synchronization modes.

One-dimensional crystals of passively-driven particles in microfluidic channels exhibit collective vibrational modes reminiscent of acoustic ‘phonons’. These phonons are induced by the long-range hydrodynamic interactions among the particles and are neutrally stable at the linear level. Here, we analyze the effect of particle activity – self-propulsion – on the emergence and stability of these phonons. We show that the direction of wave propagation in active crystals is sensitive to the intensity of the background flow. We also show that activity couples, at the linear level, transverse waves to the particles' rotational motion, inducing a new mode of instability that persists in the limit of large background flow, or, equivalently, vanishingly small activity. We then report a new phenomenon of phonons switching back and forth between two adjacent crystals in both passively-driven and active systems, similar in nature to the wave switching observed in quantum mechanics, optical communication, and density stratified fluids. These findings could have implications for the design of commercial microfluidic systems and the self-assembly of passive and active micro-particles into one-dimensional structures.

An important but missing component in the application of the kernel independent fast multipole method (KIFMM) is the capability for flexibly and efficiently imposing singly, doubly, and triply periodic boundary conditions. In most popular packages such periodicities are imposed with the hierarchical repetition of periodic boxes, which may give an incorrect answer due to the conditional convergence of some kernel sums. Here we present an efficient method to properly impose periodic boundary conditions using a near-far splitting scheme. The near-field contribution is directly calculated with the KIFMM method, while the far-field contribution is calculated with a multipole-to-local (M2L) operator which is independent of the source and target point distribution. The M2L operator is constructed with the far-field portion of the kernel function to generate the far-field contribution with the downward equivalent source points in KIFMM. This method guarantees the sum of the near-field & far-field converge pointwise to results satisfying periodicity and compatibility conditions. The computational cost of the far-field calculation observes the same O(N) complexity as FMM and is designed to be small by reusing the data computed by KIFMM for the near-field. The far-field calculations require no additional control parameters, and observes the same theoretical error bound as KIFMM. We present accuracy and timing test results for the Laplace kernel in singly periodic domains and the Stokes velocity kernel in doubly and triply periodic domains.

The number of microbial and metagenomic studies has increased drastically due to advancements in next-generation sequencing-based measurement techniques. Statistical analysis and the validity of conclusions drawn from (time series) 16S rRNA and other metagenomic sequencing data is hampered by the presence of significant amount of noise and missing data (sampling zeros). Accounting uncertainty in microbiome data is often challenging due to the difficulty of obtaining biological replicates. Additionally, the compositional nature of current amplicon and metagenomic data differs from many other biological data types adding another challenge to the data analysis. To address these challenges in human microbiome research, we introduce a novel probabilistic approach to explicitly model overdispersion and sampling zeros by considering the temporal correlation between nearby time points using Gaussian Processes. The proposed Temporal Gaussian Process Model for Compositional Data Analysis (TGP-CODA) shows superior modeling performance compared to commonly used Dirichlet-multinomial, multinomial, and non-parametric regression models on real and synthetic data. We demonstrate that the nonreplicative nature of human gut microbiota studies can be partially overcome by our method with proper experimental design of dense temporal sampling. We also show that different modeling approaches have a strong impact on ecological interpretation of the data, such as stationarity, persistence, and environmental noise models. A Stan implementation of the proposed method is available under MIT license at

We introduce methods for large-scale Brownian Dynamics (BD) simulation of many rigid particles
of arbitrary shape suspended in a fluctuating fluid. Our method adds Brownian motion to the rigid
multiblob method [F. Balboa Usabiaga et al., Commun. Appl. Math. Comput. Sci. 11(2), 217-296
(2016)] at a cost comparable to the cost of deterministic simulations. We demonstrate that we can
efficiently generate deterministic and random displacements for many particles using preconditioned
Krylov iterative methods, if kernel methods to efficiently compute the action of the Rotne-Prager-
Yamakawa (RPY) mobility matrix and its “square” root are available for the given boundary conditions.
These kernel operations can be computed with near linear scaling for periodic domains using the
positively split Ewald method. Here we study particles partially confined by gravity above a no-
slip bottom wall using a graphical processing unit implementation of the mobility matrix-vector
product, combined with a preconditioned Lanczos iteration for generating Brownian displacements.
We address a major challenge in large-scale BD simulations, capturing the stochastic drift term that
arises because of the configuration-dependent mobility. Unlike the widely used Fixman midpoint
scheme, our methods utilize random finite differences and do not require the solution of resistance
problems or the computation of the action of the inverse square root of the RPY mobility matrix. We
construct two temporal schemes which are viable for large-scale simulations, an Euler-Maruyama
traction scheme and a trapezoidal slip scheme, which minimize the number of mobility problems to
be solved per time step while capturing the required stochastic drift terms. We validate and compare
these schemes numerically by modeling suspensions of boomerang-shaped particles sedimented near
a bottom wall. Using the trapezoidal scheme, we investigate the steady-state active motion in dense
suspensions of confined microrollers, whose height above the wall is set by a combination of thermal
noise and active flows. We find the existence of two populations of active particles, slower ones closer
to the bottom and faster ones above them, and demonstrate that our method provides quantitative
accuracy even with relatively coarse resolutions of the particle geometry.

The cellular cytoskeleton is an active material, driven out of equilibrium by molecular motor proteins. It is not understood how the collective behaviors of cytoskeletal networks emerge from the properties of the network's constituent motor proteins and filaments. Here we present experimental results on networks of stabilized microtubules in Xenopus oocyte extracts, which undergo spontaneous bulk contraction driven by the motor protein dynein, and investigate the effects of varying the initial microtubule density and length distribution. We find that networks contract to a similar final density, irrespective of the length of microtubules or their initial density, but that the contraction timescale varies with the average microtubule length. To gain insight into why this microscopic property influences the macroscopic network contraction time, we developed simulations where microtubules and motors are explicitly represented. The simulations qualitatively recapitulate the variation of contraction timescale with microtubule length, and allowed stress contributions from different sources to be estimated and decoupled.

A body submerged in active matter feels the swim pressure through a kinetic accumulation boundary layer on its surface. The boundary layer results from a balance between translational diffusion and advective swimming and occurs on the microscopic length scale $$\lambda^{-1} = \delta/\sqrt{2[1 + \frac{1}{6}(\ell/\delta)^2]}$$. Here $$\delta = \sqrt{D_T\tau_R}$$, $$D_T$$ is the Brownian translational diffusivity, $$\tau_R$$ is the reorientation time and $$\ell = U_0\tau_R$$ is the swimmer's run length, with $$U_0$$ the swim speed. In this work we analyze the swim pressure on arbitrary shaped bodies by including the effect of local shape curvature in the kinetic boundary layer. When $$\delta\ll L$$ and $$\ell \ll L$$, where $$L$$ is the body size, the leading order effects of curvature on the swim pressure are found analytically to scale as $$J_S\lambda\delta^2/L$$, where $$J_S$$ is twice the (non-dimensional) mean curvature. Particle-tracking simulations and direct solutions to the Smoluchowski equation governing the probability distribution of the active particles show that $\lambda\delta^2/L$ is a universal scaling parameter not limited to the regime $$\delta, \ell\ll L$$. The net force exerted on the body by the swimmers is found to scale as $$\bF^{net} /\left(n^\infty k_sT_s L^2\right) = f(\lambda\delta^2/L)$$, where $$f(x)$$ is a dimensionless function that is quadratic when $$x\ll1$$ and linear when $$x\sim 1$$. Here, $$k_sT_s = \zeta U_0^2\tau_R/6$$ defines the `activity' of the swimmers, with $$\zeta$$ the drag coefficient, and $$n^\infty$$ is the uniform number density of swimmers far from the body. We discuss the connection of this boundary layer to continuum mechanical descriptions of active matter and briefly present how to include hydrodynamics into this purely kinetic study.

Chemically active Brownian particles with surface catalytic reactions may repel each other due to diffusiophoretic interactions in the reaction and product concentration fields. The system behavior can be described by a “chemical” coupling parameter $$\Gamma_c$$ that compares the strength of diffusiophoretic repulsion to Brownian motion, and by a mapping to the classical electrostatic one component plasma (OCP) system. When confined to a constant-volume domain, body-centered cubic (bcc) crystals spontaneously form from random initial configurations when the repulsion is strong enough to overcome Brownian motion. Face-centered cubic (fcc) crystals may also be stable. The “melting point” of the “liquid-to-crystal transition” occurs at $$\Gamma_c \approx 140$$ for both bcc and fcc lattices.