The method of MPCD was originally proposed by Malevanets and Kapral and termed stochastic rotation dynamics (SRD).^[12,13] It solves the Navier–Stokes equations on a coarse-grained level and naturally includes thermal fluctuations.^[10,11] This is achieved by modeling the fluid by N_p point-like, effective fluid particles at positions r and velocities v and each of mass m₀ and kept at temperature T₀.

At the beginning of the simulation the fluid particles are randomly placed in the simulation box and their velocities are drawn from a Gaussian distribution with standard deviation with k_B the Boltzmann constant. Then they perform alternating streaming and collision steps which are sufficient to solve the Navier–Stokes equations on a coarse-grained level because of momentum conservation in the collision step.

Figure 1 shows a two-dimensional (2D) sketch of the streaming and collision step. In the streaming step the particles move ballistically for a time Δt and the positions are updated to (Fig. 1(a))

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. The 2D sketch of the MPCD method (from Ref. [6]). (a) In the streaming step the fluid particles move ballistically. (b) In the collision step the fluid particles are sorted into cubic cells and exchange momentum with all other particles in the cell while keeping the local total momentum fixed.

Typically fluid particles and immersed objects such as microswimmers move with different time steps. While the fluid particles change their positions and velocities in relatively large steps Δt, the immersed objects update their positions – and conformation for flexible objects – at smaller time steps Δt in order to resolve the molecular dynamics (MD). In these hybrid MPCD–MD simulations in general there exist two ways to couple the immersed objects to the MPCD fluid.

First, large solid objects are usually not allowed to be penetrated by fluid particles. Momentum and angular momentum between the objects and the fluid are exchanged in the streaming step. When a fluid particle hits an object it is bounced back and the local momentum difference is transferred to the object and also modifies the angular momentum of the object.^[14] This procedure ensures translation and rotation of the objects. In addition, in order to ensure a no-slip boundary condition momentum and angular momentum can also be exchanged in the collision step by using so-called virtual fluid particles placed in a layer inside the objects (Fig. 2(a)). These virtual particles are randomly placed at each time step and participate in the collision step. The local momentum transfer to the virtual particles then contributes to the change of momentum and angular momentum of the object.

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. Coupling between the MPCD fluid and embedded objects. (a) 2D sketch of collision step. Here immersed solid spheres are shown in red and a wall in grey. Small dots in the white region indicate the fluid particles and dots in the spheres and in the wall are virtual particles, which are located in a thin layer of thickness dW = a0 in the wall and in a thin shell of thickness in the spheres. The dashed lines show the collision cell grid (from Ref. [14]). (b) 2D sketch of model polymers in MPCD fluid: polymers (light and dark green), MPCD fluid particles (black dots). The rectangular collision cell grid is shifted randomly in every time step.

Second, slender and flexible objects such as flagella, cilia, and polymers can be coupled to the fluid solely in the collision step. This can be achieved when a filament – or a membrane – is discretized by spheres connected by springs, and may be connected by other two or more body forces such as bending or twist. Figure 2(b) shows a 2D sketch of two short filaments suspended in an MPCD fluid. While the monomers have a certain size σ in order to prevent overlap of the filaments, the fluid particles are allowed to penetrate the monomers. With this respect the monomers are treated as point-like particles, which exchange momentum with the fluid particles only in the collision step. This is achieved by assigning a constant mass m_p to each of the monomers, which is typically equivalent to the mass of fluid, and a monomer occupies one collision cell.^[10] In this sense fluctuations and hydrodynamic interactions between the monomers are included efficiently. While this method was first used to study polymers in solution,^[15] it can also be applied to filaments of flexible microswimmers.

Using the method of MPCD to model the motion of microswimmers has first been applied to the dynamics of flexible swimmers such as sperm cells which consist of a spherical head, a hook, and a flexible, actively deforming tail (Fig. 3(a)). Elgeti and coworkers coupled a periodically actuated model sperm to the MPCD fluid in order to study the 3D motion of swimming sperm cells in a background fluid.^[16,17] An even simpler model in 2D has been used to study their cooperative motion, and showed the importance of hydrodynamic interactions in sperm synchronization and clustering,^[18] and the importance of hydrodynamic interaction between sperm cells and nearby solid surfaces.^[19] Recent work demonstrated how sperm can be sorted in microfluidic chambers including an array of periodic pillars.^[20] Model sperm swimming in modulated channels can be guided by optimizing the channel geometry.^[21]

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. Sketch of different microswimmers studied with MPCD. (a) Model sperm cell consisting of a solid head, an attached flexible hook, and a beating elastic tail. (b) Beating elastic filament. (c) Model bacterium consisting of a solid head with attached helical flagella driven by a motor torque connecting cell body and flagellum. (d) Active rod driven by tangential forces which are counterbalanced by forces on the fluid. (e) Slip-surface driven microswimmers such as squirmers or phoretically driven artificial microswimmers.

It is known that flexible waving sheets and filaments can swim force-free without an attached head^[4,5,22] (Fig. 3(b)). In 2D the simplest implementation is to apply a bending wave on a one-dimensional filament.^[23,24] The large scale collective properties of such waving filaments show cluster formation and collective swarming.^[23] When breaking the left–right symmetry of these filaments they are able to swim in circles, similar as observed experimentally for swimming sperm cells near solid surfaces. The hydrodynamic interactions then determine the formation of vortex arrays in collectively swimming filaments,^[25] in accordance with experimental observations.^[26] The dynamics of a single filament in confinement and periodic pillar arrays show a large speed enhancement at certain obstacle configurations, where both steric and hydrodynamic interactions are included.^[24]

Microswimmers such as Paramecia swim by synchronously beating thousands of cilia attached to their cell body. A large-scale implementation of periodic arrays of cilia shows that hydrodynamic interactions between nearby cilia determine their synchronization and the formation of metachronal waves^[27] which enable ciliates to swim.

Other implementations of flexible microswimmers in MPCD are based on classical simple model swimmers. First, it is known that a rigid one-hinge swimmer consisting of two rigid arms cannot swim at low Reynolds number because it can only perform time-reversal, reciprocal motion.^[3] This time-reversibility can be broken by using finite Reynolds numbers, non-Newtonian fluids,^[28] or flexible arms, as recently demonstrated with 2D MPCD simulations.^[29] Second, a very simple model microswimmer which is able to swim is the so-called three-sphere swimmer which periodically deforms the length of two arms connecting three aligned spheres in a non-reciprocal manner.^[30] This model swimmer has been implemented in the MPCD fluid and extended to simulate the motion of N-sphere swimmers, and to the motion of spheres which are not aligned.^[31]

Many bacteria are able to swim by rotating helical flagella attached to the cell body^[32] (Fig. 3(c)). The hydrodynamic interaction between several flagella and their flexibility allows bundle formation and synchronous rotation. In MPCD the hydrodynamic synchronization of anchored rotating flagella has been studied by discretizing a helix using beads connected by springs and including a bending and torsion potential, in order to include flexibility. These model flagella are then driven by a torque at the base, similar as for swimming bacteria, and hydrodynamic flagella synchronization can be studied by coupling to the MPCD fluid.^[33,34] The flow field of an anchored bacterium consisting of a fixed head and an attached rotating flagellum can be determined with MPCD, and it can be shown that the resulting force monopole flow field is mainly determined by the rotating flagellum, and only slightly modified by the presence of the head.^[35]

Hu et al. constructed a more refined mechanical model of bacterial flagella in order to study the bundle formation and swimming of multiflagellated E. coli bacteria.^[36] They were able to determine the influence of the specific flagella bundle on the bacterium speed, and to precisely measure the flow field. It is known that swimming E. coli bacteria close to solid surfaces swim in circles due to the counterrotation of head and flagella tail.^[37,38] Interestingly, using MPCD it can be shown that by gradually increasing the surface slip the shape of the circles is modified and eventually the sense of rotation is inverted^[39] as known for stress-free surfaces.^[40]

We have recently studied the dynamics of swimming model bacteria in explicitly modeled polymer solutions by coupling both the dynamics of the bacterium and the bead-spring polymers to the background MPCD fluid^[41] (Fig. 4). This allowed us to precisely measure the distribution and conformation of the polymers around the bacterium, and to measure the flow fields for a broad range of polymer densities and polymer types. The data from the MPCD simulations can then be used to show that a non-uniform distribution of polymers around the bacterium, and in particular around the flagella, leads to local viscosity gradients. This allowed us to demonstrate that an effective enhanced translation–rotation coupling counterintuitively increases the speed of the bacterium, although the viscosity of the fluids increases with the addition of polymers.^[41] In contrast, passive driven spheres, rods, and ellipsoids are never faster than in water, which we demonstrated recently using MPCD.^[42]

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Sketch of a model bacterium swimming in a polymer solution consisting of bead-polymers and the background MPCD fluid (not shown).

The most complex microswimmers so far implemented with the MPCD method are so-called trypanosomes which are highly deformable pathogens which cause the sleeping sickness. They are able to deform their cell body by the periodic deformation of a firmly attached flagellum, which leads to nonreciprocal cell body deformations, that allows them to swim.^[43–46]

First the dependence of the swimming speed on the mechanical properties of the organism was demonstrated,^[43–45] and later a refined model was used to quantitatively match the complex swimming patterns of experimentally studied swimming African trypanosomes.^[46]

The squirmer is a simple model microswimmer which originally has been used to model the locomotion and hydrodynamics of swimming ciliated microorganisms such as Volvox or Paramecia. Cilia at the cell surface of biological microswimmers induce an effective slip velocity close to their cell surface which propels the cell forward.^[47] In the simplest squirmer model this slip velocity is static and axisymmetric at the surface of a spherical particle (Fig. 3(e)), and analytic solutions for the flow fields and the swimming speed exist.^[48,49] By tuning a single parameter, pushers, pullers, and neutral squirmers can be simulated. Recently the squirmer became a very successful model for understanding the hydrodynamic effects in spherical microswimmer suspensions. The details of the squirmer model and the implementation with MPCD can be found in Ref. [14].

The squirmer has first been coupled to the MPCD fluid by Downton and Stark^[50] who were able to approximately recover the analytic bulk flow fields, and to determine the influence of Brownian motion on the dynamics of the microswimmer, which performs a persistent random walk. Later the rotational noise has been captured very accurately by using an angular-momentum conserving version of the MPCD method^[51] in order to study the stochastic pairwise interaction of two nearby swimming squirmers.

The squirmer can also be used to study the hydrodynamic interactions of spherical microswimmers with bounding surfaces.^[52] It can be shown that the interplay between near-field hydrodynamic interactions and thermal noise determines the time microswimmers stay at surfaces after collision. When a squirmer is confined between two parallel plates in a narrow slit, pushers, pullers, and neutral swimmers move very differently, again due to hydrodynamic near-field interactions,^[14] which can be explained by an analytic model in the noise-free limit (Supplementary Information of Ref. [53]). While the pairwise swimming of squirmers in bulk is unstable, bounding surfaces can stabilize the pair-formation of spheroidal squirmers.^[54]

Studying the dynamics of microswimmers such as bacteria in microfluidic devices enables to study active transport processes under external flow conditions. Interestingly, when hydrodynamic interactions between microswimmers and bounding walls are neglected, microswimmers do not cluster in a microfluidic channel under flow but they are homogeneously distributed. This changes for pushers and pullers which hydrodynamically interact with the channel walls, as can be demonstrated with MPCD simulations. Squirmers either cluster in the center of the channel or near the walls (pullers), or perform stable van-der-Pol-like oscillations in the channel.^[55]

Understanding the collective motion of microswimmers has become an important topic in the field of active matter physics.^[6,56] Typical model systems which neglect hydrodynamic interactions are self-propelled Brownian spheres or rods. Hydrodynamically interacting squirmers are a simple model system to study the importance of hydrodynamic interactions in collectively moving microswimmers. We studied the dynamics of multiple squirmers confined in a quasi-2D environment using MPCD^[53] and found that the type of the created flow field plays an important role in the formation of clusters and in motility-induced phase separation. Later we performed large-scale simulations in order to identify the role of the hydrodynamic pressure in the phase separation process.^[57] The influence of fluid compressibility and squirmer elongation on the collective motion in a monolayer has been demonstrated in Ref. [58].

Microswimmers which are not neutrally buoyant but, for example, heavier than water sink due to gravity, which induces external forces and Stokeslet flow fields. The competition between the sinking of heavy squirmers and the proximity of a nearby wall has been studied recently using MPCD which uncovered complex hoovering motion near walls.^[59] The collective sedimentation of squirmers under gravity shows strong convective currents and non-trivial sedimentation profiles.^[60] In the case of sufficiently low density and strong sedimentation, a monolayer of sedimented squirmers can be formed at a surface, which shows a variety of dynamic states, including a hydrodynamic Wigner fluid.^[61] When the squirmers are not only heavy, but in addition bottom heavy due to an asymmetric mass distribution, the emergent cooperative behavior shows the formation of plums, convection roles, and inverted sedimentation.^[62]

A recent study showed that the presence of polymers can enhance the rotational diffusion of squirmers,^[63] similar as observed experimentally for a spherical artificial microswimmer moving in viscoelastic fluids.^[64]

A simple model of a self-propelled rod which generates hydrodynamic flows around it has been implemented in MPCD in order to study its dynamics close to surfaces.^[65] These rods are driven by a tangential propulsion force (Fig. 3(d)) which is compensated by a force on the fluid, to ensure a force-free system. Recently the dynamics and distribution of active rods in wall-bounded simple shear flow have been analyzed with a similar model.^[66]

Schwarzendahl and Mazza implemented asymmetric microswimmers consisting of rigidly connected large and small spheres, which are able to move by imposing force dipoles on the surrounding fluid.^[67–69] They were able to simulate both pusher and puller microswimmers and the implications on their collective motion. Active dumbbells and the resulting force dipole fields have also been included in order to study the enhanced diffusion of a tracer particle in this active bath.^[70]

Artificial rigid microswimmers are typically driven by self-phoretic mechanisms such as self-diffusiophoresis, self-thermophoresis, or self-electrophoresis.^[6] These mechanisms lead to the creation of local field gradients and fluid flows tangential to the surface, which drive the particle forward.

Reactive multiparticle collision dynamics (R-MPCD) is a simulation method which is able to simulate the backcoupling between chemical concentration gradients and fluid flow using an explicit particle-based description.^[71] It combines the usual MPCD simulation technique with the ability to perform chemical reactions. Typical simulated microswimmer geometries are either Janus spheres^[25,72–76] or Janus sphere dimers.^[25,77–84] The leading order flow field of an active sphere-dimer is a force dipole, as shown in simulations^[85,86] and by an analytic calculation.^[86] This can be extended to the case where the catalytically active and inert sphere are not connected but free to move together^[87] which can lead to the active self-assembly of a microswimmer consisting of two touching spheres.^[88] Even very small Ångstrom-size motors are able to self-propel, despite of large Brownian noise.^[89,90] The formation of transient clusters appears in collectively moving chemically active Janus spheres.^[73] Cluster formation can be tuned and enhanced by modifying the particle density, and the size of the catalytic cap.^[74] More elongated collectively moving active sphere dimers in a monolayer show the formation of persistent clusters.^[91] When active sphere dimers are pinned on a 2D layer but free to rotate, strong orientational correlations exist, mainly induced by chemical interactions.^[92] If the background medium is chemically inhomogeneous, chemically active microswimmers are able to react to local gradients, or even follow chemical patterns.^[93] Time-periodic oscillations of the chemical background medium lead to a periodic dispersion-aggregation transition of active nano-dimers.^[94] The combined influence of chemical patterns and external flow can lead to a chemically induced drift in the channel.^[95] Interestingly, when the spherical symmetry is broken for active Janus particles, active rotational motion can be induced.^[96,97] When several inert and chemically active spheres are connected by flexible bonds, the explicit simulation of these active filaments shows self-propulsion, spontaneous bending, and complex shape deformations.^[98,99] By using a simple viscoelastic version of the MPCD fluid it could be shown that self-propelling dimers can enhance their speed substantially.^[100]

Self-thermophoretic microswimmers are able to self-generate local temperature gradients around their body which drive them forward,^{[84,85,101,102]} or react in thermal gradients and show phototaxis.^[103] The collective motion of chemically and hydrodynamically interacting thermophoretic sphere dimers shows collective swarming behavior.^[101]