Mesoscopic simulations at the physics-chemistry-biology interface

This review discusses the lattice Boltzmann–particle dynamics (LBPD) multiscale paradigm for the simulation of complex states of flowing matter at the interface between physics, chemistry, and biology. In particular, current large-scale LBPD simulations of biopolymer translocation across cellular membranes, molecular transport in ion channels, and amyloid aggregation in cells are described. Prospects are provided for future LBPD explorations in the direction of cellular organization, the direct simulation of full biological organelles, all the way up to physiological scales of potential relevance to future precision-medicine applications, such as the accurate description of homeostatic processes. It is argued that, with the advent of Exascale computing, the mesoscale physics approach advocated in this review may come to age in the next decade and open up new exciting perspectives for physics-based computational medicine.

Figure 1

Scales related to biological and medical applications, whose stretch can be covered by the LBPD approach. While the particle description enables the representation of complex macromolecules or cellular organizations by tracking the fate of each individual particle, the LB description is based on the collective motion of solvent molecules. The boundary between LB and PD is blurred and depends on the degree of microscopic detail required by each single application (see also Fig. 8).

Figure 2

Translocation of a biopolymer, as in the case of DNA, through a narrow pore. Different representations can be used to study multiple levels of detail during the translocation process. (a), (b) The biopolymer is represented as a simple necklace of beads by neglecting correlations stemming from the local molecular rigidity or backbone charge. (c), (d) At the next level, the macromolecule is charged and moves in an electrolyte solution, whereby a neutral solvent and counterions and coions migrate due to an externally applied electric field, giving rise to an electro-osmotic flow that ultimately causes the molecule to translocate. From [52].

Figure 3

Representation of the crowded interior of the cell as obtained from simulations. The inset illustrates the embedding of a protein in the LB mesh and each protein atom is represented via the DPM particle-fluid exchange scheme. From [16].

Figure 4

Left: Snapshot of a configuration of 1000 amyloid peptides $A{\beta }_{16-22}$ simulated in a cubic box of size $30\times 30\times 30{\mathrm{nm}}^3$ and complex flow structure generated by their motion in the surrounding solvent. Right: Evolution of the size of the peptidic aggregates as a function of time. From [160], and [43].

Figure 5

Snapshot of molecules and particles in a multiphase flow, either (a) dragging the fluid from one phase to another or (b) sitting at the fluid-fluid interface. From [190]. (c) In absence or (d) in the presence of colloids modifying the shape of the interface to a corrugated one to reflect molecular correlations. From [218].

Figure 6

Examples of standard 2D and 3D LB lattices, with 9, 19, and 27 discrete speeds, typically denoted as D2Q9, D3Q19, and D3Q27, respectively.

Figure 7

The effect of wall boundaries on the LB populations according to the bounce-back scheme. Populations hitting a wall node are reflected from the incoming fluid direction, giving rise to no-slip flow conditions.

Figure 8

The description of physical-chemical-biological (PCB) systems requires a suitable integration of microscopic details (individuality) within the universal harness dictated by symmetry principles and ensuing conservation laws which govern the macroscopic behavior (universality). Kinetic theory is expected to strike an optimal (problem-dependent) balance between universality and individuality.

Figure 9

A sketch of a high-density–low-density diffuse interface of width $w$, as typical for many biosystems (dotted line), for different values of the mean-free path, of length $\ell$. Each solid line shows a different profile of the interface, depending on its intrinsic structure. Here we assume that the wavelength of the density profile coincides with $\ell$. In the limit $\ell /w\ll 1$ (solid line, bottom), nonequilibrium effects are negligible and hydrodynamics holds. On the other hand, when $\ell$ becomes comparable to the interface width $w$ (solid line, top), nonequilibrium effects can no longer be neglected and higher-order kinetic moments must be accounted for.

Figure 10

A sketch of the fluctuating distribution $\tilde{f}$ (wiggly line), featuring rapid oscillations with respect to the nonfluctuating one $f$ (smooth line). Fluctuating hydrodynamics emerges as a consequence of the stochastic component of the distribution function.

Figure 11

Sketch of the Shan-Chen piecewise linear potential extending over the first Brillouin lattice cell at a distance $\mathrm{\Delta }x$. The potential mimics the attractive tail of a van der Waals potential (dashed line) while the repulsive one is quenched to zero. Notwithstanding the absence of hard-core repulsion, the Shan-Chen potential does not cause any unstable density pileup because, in the high-density limit $\rho \gg 1$, the generalized density $\psi (\rho )=1-e^{-\rho }$ flattops to a constant 1, yielding a zero gradient, hence zero force.

Figure 12

Ladd’s method for a suspended (spherical) body: (a) the body-fluid interface is tracked by searching the mesh points inside and outside the particle ([116]). Here a modified bounceback scheme is applied to account for the local exchange of momentum. The midpoint of the connecting links identifies the staircased surface that corresponds to the no-slip condition. (b) the particle center moves in the continuum according to the momentum exchange, and the rotational motion is accounted for by computing the total torque.

Figure 13

Examples of immersed bodies in the LBPD approach. (a)–(c) The different coupling methods between particles and fluid as detailed in the text: (a) The PPM with pointlike particles exchanging information with the nearest grid node. (b) The DPM for diffused isotropic particles. (c) The DPM for diffused anisotropic particles. (d) The time-stepping method between particle dynamics (PD) and fluid dynamics (LB). The time-stepping scheme of the LB and PD individual components can be either synchronous or asynchronous. Typically, the PD component ticks more frequently than the LB component, allowing for more efficient simulations. (e) A typical configuration of the fluid around a moving particle. Hydrodynamic response manifests itself as Stokes flow field (upper part) and long-range flow structure around a spherical particle or its far-field flow pattern (lower part).

Figure 14

Proteins immersed in a LB solvent are treated as mesoscopic objects. Consequently, their representation is coarse grained, that is, a level of detail based on a handful of pseudoatoms per amino acid (five pseudoparticles in the current case from the OPEP model) ([35]). The pseudoparticles are needed to encode the details of the peptidic backbone and effective interparticle interactions.

Figure 15

Three-level representation of blood flow. Upper panel: continuum macroscale (cm). Bottom left panel: the granular nature of red blood cells starts to be apparent (100 μm) and is fully revealed at the level of the cell spacing (10 μm, bottom right panel). From [17].

Figure 16

The road map of the muphy code, based on the Gordon Bell performance obtained over the years. From [17].

Figure 17

The flow structure around a macromolecule resolved by a multiresolution mesh being finer in proximity and inside the macromolecule.

Figure 18

Diagram for the multiscale approach where Eulerian (E) and Lagrangian (L) agents can move across scales and across representations. The boxes of different size in the left column indicate the coarse or fine grain of the Eulerian representation, and the ball-and-stick representations in the right column indicate particle-based molecules at different resolution. The arrows indicate the exchange of agents between E-E scales, L-L scales, and E-L representations.

Figure 19

Snapshot of globular proteins in solution, shown in different visual representations. (a) Proteins are shown with wirelike bundles and the surrounding hydrodynamic signal is represented via an isosurface of its velocity field. Streamlines generated by selected proteins in the crowded macromolecular environment show how the hydrodynamic disturbance propagates in the aqueous solution. (b) Proteins represented via conventional ribbon-sticks style and the constant velocity isosurface that illustrates the complexity of the flow structure.

Figure 20

Multiscale representation of vesicles transporting proteins. Vesicle firing from a membrane where the vesicles are modeled as immiscible fluid phase separating from the aqueous host. The central panel zooms the vesicle bilayer with liquid water in the interior. The right panel depicts a molecular representation of the vesicle containers with an heterogeneous protein suspension in the interior and membrane proteins embedded in the bilayer.

Figure 21

The major phases of hemostasis following the vessel injury: (a) a platelet plug is formed to rapidly stop the initial bleeding; (b) a mesh of fibrin is made at the injury point to make the platelet plug stronger; and (c) finally the clot is formed by the coagulation of fibrin, red and white blood cells.

Figure 22

CPU vs GPU optimal data layout. Row major ordering is assumed for storing multidimensional arrays in memory (typical of $C$ and C++ languages). AoS stands for array of structures; SoA stands for structure of arrays.

Figure 23

Communication schemes: (a) with no overlap between communication and computation and (b) multi-GPU scheme using two streams.

Figure 24

Organization of mesh nodes according to the indirect addressing scheme for an irregular domain. (a) In the D2Q9 scheme, red (dark gray) squares correspond to fluid regions and green (light gray) squares to wall regions. Inactive nodes are in white and are not numbered since they are not stored in memory. Entries in the connectivity matrix are shown for fluid node 9. (b) Indirect addressing on the GPU for a 3D case. Each GPU thread handles a subset of fluid nodes and a subset of populations.