Meso-scale Algorithm for Flowing Nematics

Most of my projects investigate some microscopic object (or objects) such as polymers or microbial swimmers immersed in a simple background fluid. Since it’s the dynamics of the objects that we’re typically interested in, we treat the fluid as simply as possible. In order to do this, we often use the MPCD method (described here).

However, not all fluids are simple. Some biological fluids have intrinsic elasticity, while others are made of or include rod-like molecules. These fluids are called “nematic” — from the Greek word for “threadlike”. Since the molecules are rod-like, they can align; and when they align they are partially ordered. Ordered solids are called crystals and so these fluids are called liquid-crystals. Since many biological molecules are rod-like, such as filamentous-actin, their solutions can be liquid crystalline.


Dealing with flowing liquid crystals is notoriously difficult. When the liquid crystal is just the background fluid to the dynamics of the colloids, polymers, or microbes, significant computational effort must be spent to simulate the background. This is why we developed a computationally efficient and versatile multi-particle collision dynamics (MPCD)-based algorithm to simulate nematohydrodynamic baths.


Our key realization was that the MPCD algorithm isn’t limited to solving the Navier-Stokes equation for momentum transport in a fluid, but can really solve any transport equation so long as a suitable “collision” operation (now more generally referred to as an “exchange”) can be defined and the continuum field of interest can be discretized as a particle property. We gave each point-like MPCD particle an orientation that would interact by an order-collision operation based on Maier-Saupe theory for liquid crystals in equilibrium.

This simple algorithm works very well. It is accurate in 2D and 3D; it reproduces defect dynamics; it has an intrinsic elasticity; it can have tumbling or shear aligning regimes when it flows; and it can handle complex boundary conditions and mobile macromolecules, colloids and swimmers.

Lattice Boltzmann Simulations of Swimming of Microbial Organisms

Microorganisms live in a very different world than us — a world commonly submerged in water. But, at such microscopic sizes, swimming in water seems far more like swimming through tar. Understanding the dynamics of how microorganisms travel through liquids can help us control their spread.

Since the microbial world is a liquid one, flows play an important role in setting microbe motion. Bacteria, for instance, commonly have “run-tumble” dynamics: They swim forward for a time then briefly tumble in a new, random direction. Walls and obstacles can impact these dynamics, not only by simply getting in the way but also by changing the flows caused by the bacteria. In this way, even distant objects interact with the swimmers through flow.

Images by Joost de Graaf.

These hydrodynamic interactions are challenging to predict. We’ve worked with Joost de Graaf (University of Edinburgh), who has developed versatile models for computationally simulating the dynamics of self-propelled particles (such as living bacteria) within fluids. These simulations build swimmers out of beads, which appear raspberry-like in shape. The resulting models reproduce well the hydrodynamics of swimming cells, like E. coli or sperm, which means that in the future these models may be able to accurately predict swimming dynamics in complex geometries and flows.

Microbial Swimming in Films

Living near interfaces can be beneficial to microbes. Liquid-air interfaces may be the source of oxygen and sunlight. Solid surfaces accumulate sediments, which can include valuable nutrients. They can also offer anchoring points for the formation of protective biofilm. Films of liquid over a solid wall have both of these types of beneficial interfaces. However, an unavoidable consequence of being near surfaces is a shearing when the film flows. These flows can have a large impact on the swimming dynamics of microbes.

Image by Arnold Mathijssen.

In a series of papers with Arnold Mathijssen and Amin Doostmohammadi, we investigated the swimming dynamics of motile microbes in flowing films. We began by building an analytical model for the flows produced when microbes swim in a film. We were able to reproduce many of the different flows generated by swimmers. These include flagellated bacteria that swim by propelling fluid back and then pushing fluid forward (called “pusher-type swimmers”), and puller-type swimmers that pull fluid toward their head and drag their body forward. We also reproduced ciliated swimmers with their many hair-like flagella.

We then put our swimmer models in a flowing film and looked at the dynamics. We found there is a sweet spot for rinsing bacteria away. When the film isn’t flowing, the microbes are attracted most strongly to the solid surface and they stay near the floor until the flow is strong enough to peel microbial swimmers off the surface. That is to say, at moderate flows bacteria in a film are washed away but that’s not the end of it. If much stronger flows are applied, the flow turns the microbes back down towards the surface, where they are more likely to form biofilms.

Image by Arnold Mathijssen.

Run-tumble Dynamics vs. Hydrodynamic Interactions in Films

Swimmers like E. coli have “run-tumble” dynamics, meaning they swim forward for a time (run stage) then rapidly and randomly change to new direction (tumble stage). While they swim in the run stage, they perturb the surrounding fluid and interact with their surroundings through the flow they generate. This causes swimmers to be attracted to walls.

We asked “When does the hydrodynamic attraction win over the randomizing run-tumble dynamics in a film?” We found if we modelled the dynamics as just an effective constant noise then the swimming strategy of the microbes makes quite a difference. However, if the swimmers have a proper run stage followed by a brief and random tumble stage then the swimming strategy is not nearly as important.

We believe this will be of interest to other scientists in our field because it means that when they make simulations of swimmers they should prioritize getting the run-tumble dynamics correct before attempting to reproduce the more complicated but less important hydrodynamics of microbial swimming.

Upstream Swimming in Complex Fluids

We’ve worked on computational and theoretical models for swimming microbes and looked at their dynamics in biologically relevant environments, such as flowing films. Another important consideration is the nature of the fluid that the microbes swim through. An important aspect of many fluids is that they have a dual fluidic and elastic (viscoelastic) nature. Biological examples of such so-called non-Newtonian fluids include cervical mucus, extracellular matrix surrounding biofilms, and flowing blood. To make matters more complicated, microbes such as sperm and other swimming cells must often swim upstream.

With Arnold Mathijssen and Amin Doostmohammadi, we considered a microswimmer swimming in a channel geometry and subjected to a flow. Previous studies with this set-up were done by Zottl and Stark, who established that the mathematics describing the microbe’s swimming are related to the swinging of a pendulum. We built on their work to predict that swimming cells in biological flows naturally reorient toward the centre of channels, where they migrate upstream. This is because of the elasticity in such fluids.

The reorientation also happens because the flow turns the swimmers away from the walls, like in normal, Newtonian fluids (i.e. water). We found that swimmers in fluids that become more viscous when sheared (like cornstarch that thickens when stirred) swim upstream more quickly than in fluids that thin (like ketchup that spills suddenly after slowly creeping). The swimmers move upstream faster in shear-thickening flows because they can spend more time in slower flows near walls, where they have the chance to sneak further upstream.

Figure by Arnold Mathijssen.

Swimming speed is used to assess the fertility of human semen so we suspect that fertility estimates based on speed could be improved by performing these measurements in viscoelastic fluids. The sperm in such a flowing fluid would swim upstream along the centreline where they can have a fair race, without some cheating by sneaking up along the walls.

This work will soon appear in the journal Physical Review Letters.


Field-Flow Fractionation (FFF; introduced here) separates samples of particles in a flow by applying an external field that causes distinct solutes to have different concentration distributions. The greater the difference in the distributions, the greater the separation of different solutes.

Traditionally, FFF strives to have just a single external force applied to the solutes — more than one external field has been assumed to (a) just increase the net forces when the fields align (termed additive-mode FFF by Bruce Gale) or (b) simply decrease the net forces when the fields compete. One would never apply opposing forces on purpose.


However, while trying to develop models for understanding hyperlayer-FFF (when unavoidable lift forces oppose the external field), I was surprised to find that two opposing fields can produce a dramatic drop in elution time of very specific solute sizes. Gary Slater and I termed this novel technique adverse-mode FFF”. By applying two strongly-opposing external forces, which scale differently with solute size, a sharp peak occurs in the speed of a solute for a given size. For solutes smaller than the peak position, the force that scales less strongly with size dominates. For larger solutes, the force that scales more strongly with size dominates. Exactly at the peak position, the two opposing forces balance. By increasing both field strengths together, the peak width can be narrowed and be made quite sharp.


All of this occurs for realizable field strengths. To measure an unknown solute size, adverse-mode FFF could following this recipe:

  • An external field that scales directly with solute radius (such as flow-FFF, electrical-FFF or thermal-FFF) is applied to an eluting sample.

  • A second field that scales with the solute volume (radius cubed; such as centrifugal-FFF) is applied to oppose the first field.

  • By incrementally decreasing the strength of the second field from an initial large value, the peak position would shift to smaller values and an abrupt increase in the speed of the solute would occur when the peak location matches the solute’s size.

  • The solute size can then be calculated.

  • Having found the solute size, the accuracy can be improved by increasing both external field strengths together.

DNA in Nanoslits

DNA doesn’t just float around. It it tugged on and wrapped up by the protein machinary within cells. All of this pushing and pulling occurs within the confined space of the nucleus. In artificial systems, experimentalists pull on DNA and other semi-flexible biomolecules, stretching them within channels or slits.

Simulation snapshots by Hendrick de Haan.

However, the model relating the extension of DNA to the amount of force applied is for fully 3D or perfectly 2D, like in a film. In strongly confining slits, the walls make it easier to stretch the DNA but not as easy as in 2D. The confining walls make it easier because they obstruct the DNA and so limit the number of conformational states that the DNA can have. Since the number of states sets the entropy and the entropy works against the applied stretching force, confinement makes stretching easier.

With Hendrick de Haan (University of Ontario Institute of Technology), we developed a theory for describing the extension of DNA in a channel of any height as a function of force. This theory is built on the classical theory of Marco-Siggia for DNA in 2D and 3D and uses an effective dimensionality between 2 or 3 when the forces are not too strong. When the forces are very strong a mathematical approximation must be applied for strong and weak confinements.

Entropic Collapse of Coarse-Grained Chromosomes

Figure 1
Illustration by David Goodsell from Chromosome, Cell Cycle, and Entropy.

Entropy can be tricky to understand, but one rule often quoted in pop-culture is that entropy increases as time goes on, which leads to more and more disorder. Depletion forces can therefore be surprising since they seem to break this rule-of-thumb; the depletion force is an entropically induced attraction between large colloids due many smaller surrounding “depletant” particles. By bringing the large colloids together, the attraction orders the larger colloids but this is at the benefit of the many small depletants, which gain increased entropy and are more disordered.


Such entropic forces exist in biological systems. For example, depletion forces lead to fibre bundling, aggregation of red blood cells, and change the reaction kinetics of polymerase chain reactions. They are also ubiquitous in the complex and crowded interior of cells.

Prior to our work, it was hypothesized that depletion forces acting on chromosomes in simple bacteria by surrounding cytoplasmic proteins could explain the collapse from a swollen conformation to a compact state. We performed simulations to explicitly test whether depletion-induced attraction is actually enough to collapse model chromosomes. The resulting simulations gave us strong evidence that depletion effects are enough to collapse chromosomes; however, this collapse is a continuous phase transition rather than a discontinuous one (as is seen in experiments). This suggests that some proteins must act in specific biological ways (rather than as generic depletants) to fully account for the compaction of nucleoids.


To read more, see our paper on the entropic collapse of coarse-grained chromosomes, or a News and Views by Suckjoon Jun discussing our work.


Slip-FFF with Super-Hydrophobic Surfaces

The flexibility of FFF stems from the large set of transverse external fields that can be used. Gravitational fields, centrifugation, electrical fields, magnetic fields, dielectrophoretic fields, acoustic fields, photophoretic fields, cross-flows (both symmetrical and asymmetrical) and thermal gradients have all been used to generate non-uniform solute concentration distributions.

But FFF also depends on a non-uniform flow profile. It is generally taken for granted that the profile in a channel obeys Poiseuille flow. By changing the flow profile we may change the retention ratio (the mean velocity of the solutes normalized by the mean fluid flow).

One way to alter the flow profile is to have a slip wall. Usually, fluids flowing over a solid surface obey a no-slip boundary condition, meaning that the speed goes to zero. However, modern nanoengineering can construct ultrahydrophobic surfaces by nanopatterning structures on the surface that reduce the drag on the fluid.

Another idea is to superimpose Poiseuille flow and electroosmotic flow with a thin Debye layer. This would create a non-zero velocity at the wall (called the Smoluchowski slip velocity, if you’d like to look it up).

With these motivations in mind, we investigated how the retention ratio would change for various slips. We found an number of interesting conclusions:

  1. Hydrodynamic chromatography never benefits from slip.
  2. In a moderate external field Normal-Mode FFF can be somewhat improved by slip at the depletion wall (the top wall).
  3. Slip at the depletion wall typically increases the range of Steric-Mode FFF with little impact on resolution.
  4. Normal- and Steric-Mode FFF persist for weaker fields when there is slip at the depletion wall.

But the most interesting result by far is that we once again found a novel operational-mode. If the accumulation wall (the bottom wall) has an extremely large slip and the external field is strong then the retention ratio rises surprisingly rapidly over a small range for the tiniest particles. Mathematically, this looks like the hydrodynamic chromatography regime of FFF but the particle size range and the resolution are more comparable to Normal-Mode FFF so we call this new mode slip-FFF.

Find more details or the citation for this in our article in the Journal of Chromatography A.

Polyelectrolyte Brushes

Nanoengineering surfaces is a concept that comes up at multiple occasions within these research projects. I’m fascinated by the idea of constructing something that is nanoscopic in at least one dimension that can have macroscopic effects on the operation of a technology. In all these projects, the idea is to use nanoengineered surfaces as components in lab-on-a-chip systems.

One appealing idea is the concept of smartnanochannels, i.e. channels that would respond to stimuli such as variation in temperature, pH, or concentrations of specific chemicals. The most obvious surface modification that could have the potential to respond to such signals is grafting (or otherwise associating) polymer chains to interfaces. In fact, it is charged polyelectrolyte brushes that have the greatest potential as smart brushes. Because of the long-range electrostatic coupling between chains and free counterions, the structure of polyelectrolyte brushes depends sensitively on many factors, which could be used to switch the brush from one state to another.

Swelling and deswelling by an external voltage, applied normal to the grafting surface has been experimentally demonstrated and suggests a convenient and quick switch that acts on both the polyelectrolyte brush and also the counterions suspended in the solvent. The extension/collapse transition could act as a stimuli responsive gate, limiting the passage of flowing fluid or ions and so has potential in the design of functional and controllable MEMS devices.

With our collaborators Dr.Pai-Yi Hsiao and his student Yu-Fan Ho from the National Tsing Hua University, we have found some very interesting results from studying this system. More to come once we know where our manuscript stands. For now you can read more about the electrohydrodynamics of soft surfaces in our review on electrophoresis.

Electrophoresis Algorithm with Finite Debye Layers

We have developed a hybrid MD-MPCD algorithm to simulate the electrophoresis of polyelectrolytes in a way that coarse-grains electrostatic interactions.

The Problem

A charged monomer in solution attracts oppositely charged ions out of the fluid to form a counterion sheath. This is called the diffuse Debye layer and the charge density of the Debye layer obeys the Poisson-Boltzmann equation for which there is no analytical solution for a spherical monomer. If the potential is rather small then the Debye layer can be described by the Debye-Hückel approximation (the counterion density decays exponential over the characteristic Debye length).

Beyond the Debye length, electrostatic interactions are screened by the oppositely charged cloud of counterions. Perhaps more surprisingly, hydrodynamic interactions due to electrophoretic motion are also screened. To understand why imagine that the monomer moves under the action of an electric force. It shears the fluid, pulling it along with it. But the counterion sheath has an equal but opposite charge. The field forces it in the opposite direction and it also shears the fluid. The shears (almost) cancel out. (I throw in the “almost” because there is a rapidly decaying shorter-ranged field that remains).

If you were reading carefully you noticed that I explicitly said that “hydrodynamic interactions due to electrophoretic motion are screened”. Normal, random motion from thermal noise or drag from pressure driven flows are not screened. When the polymer moves, its counterion cloud moves along with it (not against it as it did for an electric field). This means that the electrophoretic mobility looks like a free-draining (or Roussian) polymer but the diffusion coefficient (hydrodynamic mobility) looks like a hydrodynamically oblique (Zimmian) polymer.

Electrophoresis: When Hydrodynamics Matter

So you see that hydrodynamic interactions are usually screened in charged polymers. But this isn’t always the case. You can read all about cases when this screening isn’t perfect in our review if you are interested.

In order to study systems in which hydrodynamics aren’t perfectly screened a coarse-grained algorithm that can properly account for the electrohydrodynamics is needed. Usually simulations that account for both hydrodynamics and electrostatics are compuationally expensive (both hydrodynamics and electrostatics are long-range and computationally intense) or are only valid in the infinitely thin Debye layer limit. We developed an algorithm for the electrohydrodynamics of polyelectrolytes that is less computationally intense but is still accurate for finite Debye layers.

We embed a charged MD polymer in an MPCD fluid (as a detail the MD beads are included in the MPCD collision step to couple the polymer to the fluid). At each time step, we use the counterion distribution from the Debye-Hückel approximation to bestow MPCD particles with a temporary charge. When an electric field is applied the MD beads are forced in one direction (electrophoresis) and their not-explicitly-included-counterions are forced in the opposite direction (electro-osmosis). During the MPCD collision step this momentum is mixed – the MD beads shear the fluid and the electro-osmotic flow drags on the MD beads. For long polymers in free-solution this cancels out so that the mobility is length independent. For short chains, the experimentally observed rise is obtained (one point that I didn’t mention is that MPCD particles can get charge from multiple monomers and if their charge becomes to great they are said to “condense” to the chain, which is to say they loose all their charge and the monomers that they got their charge from also loose it. This Manning condensation is very important but if you are worried about these sorts of details you really should go read our paper).

Now we use this algorithm to study situations in which we are either unsure whether hydrodynamic interactions are totally screened or the impact of partial screening. This work was extended to look at electro-osmotic flow generated by tethered polyelectrolytes when I visited Dr. Pai-Yi Hsiao at the National Tsing Hua Univeristy in Taiwan. Currently, Martin Bertrand and I are using it to investigating some surprising effects of confinement on electrophoresis.

This algorithm was developed with Owen Hickey.