The structure and dynamics of confined suspensions of particles of arbitrary shape are of interest in multiple disciplines from biology to engineering. Theoretical studies are often limited by the complexity of long-range particle–particle and particle–wall forces, including many-body fluctuating hydrodynamic interactions. Here, we report a computational study on the diffusion of spherical and cylindrical particles confined in a spherical cavity. We rely on an immersed-boundary general geometry Ewald-like method to capture lubrication and long-range hydrodynamics and include appropriate non-slip conditions at the confining walls. A Chebyshev polynomial approximation is used to satisfy the fluctuation–dissipation theorem for the Brownian suspension. We explore how lubrication, long-range hydrodynamics, particle volume fraction, and shape affect the equilibrium structure and the diffusion of the particles. It is found that once the particle volume fraction is greater than 10%, the particles start to form layered aggregates that greatly influence particle dynamics. Hydrodynamic interactions strongly influence the particle diffusion by inducing spatially dependent short-time diffusion coefficients, stronger wall effects on the particle diffusion toward the walls, and a sub-diffusive regime—caused by crowding—in the long-time particle mobility. The level of asymmetry of the cylindrical particles considered here is enough to induce an orientational order in the layered structure, decreasing the diffusion rate and facilitating a transition to the crowded mobility regime at low particle concentrations. Our results offer fundamental insights into the diffusion and distribution of globular and fibrillar proteins inside cells.

Diffusion under confinement is central to multiple physical, chemical, and biological systems, including colloidal and protein suspensions, devices for particle separation, and transport through membranes. A model system to study the diffusion and structure of highly concentrated particles under confinement could offer insights into the dynamics of crowded macromolecules, such as proteins, inside cells where they typically occupy 20%–40% of the cytoplasm volume.1 

Previous studies have shown that crowding between macromolecules affects reaction rates of equilibrium reactions2–6 and hinders the diffusion of intra-cellular particles.7In vivo experiments, using fluorescence recovery after photo-bleaching (FRAP) techniques, have reported that the apparent diffusion coefficient for green fluorescent proteins (GFP) in E. coli’s cytoplasm is about 11 times lower than that in water.8–10 Even though a variety of intracellular activities, namely, metabolism, cellular homeostasis, signaling, transcription, translation, and locomotion, strongly depend on diffusion of the intracellular macromolecules, the mechanisms behind the hindered diffusion are not fully understood. A review by Skolnick discussed several factors that could lead to hindrance, including the viscosity of cytoplasm, steric effects, hydrodynamic interactions (HI), and other short-range interactions between particles.11 A significant effort that relied on Brownian simulations of 50 types of macromolecules modeled as spheres inside E. coli’s cytoplasm was able to reproduce the translational diffusion coefficient of GFP.12,13 The authors demonstrated that steric repulsions cannot explain the hindered diffusion and suggested that electrostatic and other short-range interactions are essential variables to consider. However, short-range interactions modeled as van der Waals potentials can always be tuned to match experimental results, and these authors were unable to provide conclusive remarks regarding the effects of short- and long-range HI. In other work, Brownian Dynamics (BD) simulations considering fluctuating HI have been performed in the bulk,14 confirming that HI plays a nontrivial role in the hindering of macromolecular diffusion. Unfortunately, the mobility variations induced by confinement were not considered in this work. Confinement, using a network of wall particles that are constrained by a predefined potential, was included in a subsequent study to represent a cell membrane.15 In that work, the non-slip boundary condition was not strictly satisfied, allowing the flow to penetrate the membrane and driving a tangential component of the particle average velocity. From the point of view of the role of hydrodynamic interactions on mobility, previous efforts were centered on the study of dynamics for single or many particles immersed in an unconfined viscous fluid16–24 and that of particles moving near a wall25–30 or confined in a slit,26,31–34 in cylindrical geometries,26 and in a spherical cavity.25,35–37 Recently, a theoretical study examined the behavior of a concentrated colloidal dispersion of spheres confined in a spherical cavity.37 That framework relied on a set of hydrodynamic tensors that capture far- and near-field (lubrication) hydrodynamics between particles and walls. The authors studied the structure and diffusion of the system using a Stokesian dynamics (SD) approach36 and reported that the confinement, crowding, and HI collectively lead to an anisotropic micro-structure, which then induces position-dependent and anisotropic short- and long-time dynamics. That study was, however, limited to spherical particles.

We have developed an efficient computational framework to perform BD simulations of arbitrarily shaped particles confined in any type of geometry. We use an Immersed-Boundary (IB) method to represent the suspended particles, a parallel finite element general geometry Ewald-like method (pFE-GgEm)38 to calculate the confined Green’s functions, and a Chebyshev polynomial approximation to satisfy the fluctuation–dissipation theorem. In this work, we use this methodology to study how steric repulsion, short- and long-range hydrodynamic interactions, confinement, particle volume fraction, and particle shape affect the structure and the diffusion of spherical and cylindrical finite-size particles confined in a spherical cavity. The cylinders are selected to break the three-dimensional symmetry of the particles, a feature that is common in protein structures.

The paper starts with a description of the particle model, the geometry, the methods, and the dimensionless variables to characterize the system. The BD and the IB-pFE-GgEm methods are briefly summarized. We then proceed to present and discuss the results, including the spatial and orientational ordering of the particles and the short- and long-time diffusion behavior. We finish the paper with a summary of the most important findings.

Let us consider N mono-disperse and semi-rigid particles embedded in a viscous fluid of viscosity η, confined in an spherical cavity of radius R. Under a zero Reynolds number condition, the N-body force/torque balance on the particles is

FH+FB+FC+FEV+Fext=0,
(1)

where FH is the 6N hydrodynamic force/torque vector, FB is the Brownian force/torque vector, FC is the force/torque vector containing configurational forces, FEV represents the force/torque vector due to excluded volume interactions, and Fext includes all external forces/torques.

The evolution of the suspended particles, from Eq. (1), is carried out through the grand mobility or resistance tensors that relate the hydrodynamic forces/torques to the translational and rotational velocities of the particles.39–41 SD42–44 and boundary integral methods (BIM)39,45 have been used extensively to solve this “mobility problem.” The regularized Stokeslets,46 the accelerated BIM,47 and the Immersed Boundary (IB)48,49 approaches all provide computational efficiency and simplicity, typically to improve (or avoid) the calculation of the single- and double-layer hydrodynamic potentials of the suspended particles. In particular, the IB method represents the surfaces of the suspended solids as a distribution of discrete force densities that, together with a surface force description and Stokes equations, generate the temporal evolution of the suspended particles. This is the approach that we use in this work.

The surface of each suspended particle is discretized into a set of NIB nodes that constitute a mesh, similar to boundary element methods.45 On the surface nodes, we define structural spring potentials that maintain particle shape, volume, and surface. The force balance on each of the N particles is then translated into the NIB surface nodes as follows:

fνH+fνB+fνC+fνEV=0,
(2)

for every node ν = 1, …, NIB, where fνH is the hydrodynamic force, fνB is the Brownian force, fνC is the constitutive force, and fνEV is the force from all the excluded volume interactions: particle–particle and particle–wall.

Assuming that the probability density for the nodal positions is a continuous density for the Fokker–Planck equation,50 an equivalent stochastic differential equation for the motion of the nodes is written as follows:51 

dR=MF+RDdt+2BdW,
(3)

where R=(x1,x2,,xN×NIB) denotes a 3N × NIB vector containing the spatial coordinates of the nodes, D = kBTM is the 3N × 3N diffusion tensor, kB is the Boltzmann constant, T is the absolute temperature, and M is the (3N × NIB) × (3N × NIB) mobility tensor. In addition, U = M · F contains the 3N fluctuating velocities from the hydrodynamic interactions and F is the 3N × NIB vector that contains the non-HI and non-Brownian forces on the nodes. The divergence of the diffusion tensor RD is the drift resulting from the configuration-dependent mobility of the confined particles and dW is a random vector, the components of which are obtained from a real-valued Gaussian distribution with zero mean and variance dt; dW is coupled to the diffusion tensor through the fluctuation–dissipation theorem: D = B · BT.

In IB methods, the force distributions at moving solids are discretized as distributions of regularized point-forces. The “smoothing” function for the delta function scales as the distance between the surface nodes that are used to represent the moving particles. Consequently, the structural forces on the particles, fνC, define a force density as follows:

ρIBC(x)=ν=1NIBfν C(xν)δIB(a,xxν),
(4)

where δIB(a, x) is a smoothing function that is a modified Gaussian (details available in Ref. 52). The regularization parameter ξIB in δIB is related to the characteristic length h for the node spacing on the particle surface, i.e., ξIBh−1a−1. The rational behind δIB(a, x) is to ensure that the regularized force on each node is spread over the length scale of the associated surface elements, thereby preventing fluid from penetrating the particle surface.

In this work, we consider the particles as “semi-rigid,” where each node on the particle is linked to its neighboring nodes by elastic springs with a prescribed large stiffness constant. The nodes are also connected to the particle center-of-mass by an elastic spring to conserve the desired shape. For simplicity, we assumed that each link is a linear spring where the force acting on the point ν by the point μ is given by

fνμC=krνμq0xνμrνμ.
(5)

Here, k is the spring elastic constant, q0 is the equilibrium spring size for each specific situation, xνμ = xνxμ, and rνμ = |xνμ|. For each particle, a spring network is formed (see Fig. 1), resulting in an internal nodal force that resists the deformation of the particles. In addition, at each surface node of each particle, a purely repulsive Lennard-Jones (LJ) force is added to account for particle–particle and particle–wall excluded volume interactions. This LJ force comes from two contributions: (1) interactions with nodes from a different particle and (2) interactions with the wall. It is equivalent to the negative gradient of the LJ potential, which is defined by

VLJ(r)=4kBTσr12σr6+kBT
(6)

for r ≤ 21/6σ and zero otherwise. Here, r is the Euclidean distance between nodes of two different particles or between the node and walls, and σ = 2.2a is chosen empirically to guarantee that each surface node has an excluded volume of radius a. The translational and rotational velocity of particles are calculated by integrating the velocity over the surface mesh of each particle, thereby satisfying the force and torque balance. For completeness, we include a validation of our IB approach in  Appendixes A–C: (i) fulfillment of Stokes’ law of particles under confinement, (ii) validation of the fluctuation–dissipation theorem (diffusion of particles), and (iii) consistency of the particle shape as a function of the spring constants.

FIG. 1.

Snapshot of the spherical cavity of radius R containing spherical particles with ΦHI = 10%. The spherical particles radius is rS, while the size of the cylindrical particles is determined by rC and hC. The surface of the particles is given by a collection of discrete nodes that are connected to six neighbors (“trimesh”), similar to boundary element discretizations, and with a characteristic node separation of ahξIB1. The neighboring nodes are connected with surface springs (black); each node is also connected to the particle center-of-mass (red spring). A repulsive Lennard-Jones excluded volume is included on each surface node shown schematically in the particles’ cross section by the blue circles. The characteristic size of the repulsion is given by σ = 2.2a.

FIG. 1.

Snapshot of the spherical cavity of radius R containing spherical particles with ΦHI = 10%. The spherical particles radius is rS, while the size of the cylindrical particles is determined by rC and hC. The surface of the particles is given by a collection of discrete nodes that are connected to six neighbors (“trimesh”), similar to boundary element discretizations, and with a characteristic node separation of ahξIB1. The neighboring nodes are connected with surface springs (black); each node is also connected to the particle center-of-mass (red spring). A repulsive Lennard-Jones excluded volume is included on each surface node shown schematically in the particles’ cross section by the blue circles. The characteristic size of the repulsion is given by σ = 2.2a.

Close modal

In what follows, we use dimensionless variables for all results and discussions. We use a as the characteristic length scale and the nodal spacing diffusion time, a2ζ/kBT, as the characteristic time scale. We set kBT as the scale for energy and kBT/a for the force. The friction coefficient ζ is related to the fluid viscosity η and a through Stokes’ law: ζ = 6πηa, and the nodal diffusion coefficient is D0 = kBT/ζ.

Our simulation method, denoted by IB-pFE-GgEm, is an O(N) algorithm that includes hydrodynamic interactions for confined, large-scale suspensions of finite-size particles of arbitrary shape. Details can be found in our previous work.38 Briefly, the O(N) algorithm consists of three major components: (a) pFE-GgEm routines to calculate the Green’s function (Stokeslet) for any geometry, (b) Fixman’s mid-point algorithm for time integration, and (c) the Chebyshev polynomial approximation for the fluctuation–dissipation theorem. The parallel finite element GgEm (pFE-GgEm) routines are built using open source libraries, thereby facilitating usage of our software. The routines can be downloaded at http://miccomcodes.org as part of the Continuum-Particle Simulation Suite (COPSS) from the mid-west Integrated Center for Computational Materials (MICCoM).

We consider spheres and cylinders, of equal volume, that are suspended in a spherical cavity of radius R = 15. The particles’ radius is rS = 3 (volume VHI=4/3πrS3), while the cylinders’ size is determined by rC = 2.62 and hC = 2rC (volume VHI=πrC2hC). Figure 1 shows the details of the system. According to our semi-rigid particle model, there are two ways to define the particle concentration in a cavity of volume V. One uses the hydrodynamic volume fraction, ΦHI = NVHI/V, and a second one is based on the excluded volume effective size, ΦEV = NVEV/V. Each surface node has an excluded volume of radius a, thus each spherical particle has an excluded volume of VEV=4/3π(rS+a)3, and each cylindrical particle has an excluded volume of VEV=π(rC+a)2(hC+2a). We use the hydrodynamic volume fraction as the relevant scale for particle concentration. In this work, ΦHI = [5%, 10%, 15%, 20%], which is equivalent to ΦEV = [12%, 24%, 36%, 48%]. The lower volume concentration selected here is 5%, which represents a limiting volume fraction between the dilute and finite concentration regimes. We found that all results for dilute systems, starting from infinite dilution, are close to those for the limiting volume fraction of 5%. We start our analysis by exploring the structure of the particles using the particle number density as a function of radial position within the cavity. The number density (the probability that a particle is at a specific location) is calculated by discretizing the spherical cavity into m bins (spherical shells) with an even spacing in the radial direction. The shell radius of the ith bin is bi = (i + 0.5)R/m. The particle number density n(ri) = ⟨N(ri)/Vi⟩, where N(ri) is the number of particles at shell that is located a distance ri from the center of the cavity, Vi is the volume of the shell, and ⟨ ⟩ represents an ensemble average over time. Figure 2 shows the number density of the particles within the cavity as a function of particle concentration. Figure 2 (top) is for spherical particles, while Fig. 2 (bottom) is for the cylindrical ones. At low concentrations, ΦHI = 5%, the number density is uniform throughout the cavity and goes to zero once the particles are in contact with the wall. Note that the maximum density is at r ∼ 10.5, which is smaller than wall contact. In our model, particles can never “touch” the wall because of the strong repulsive Lennard-Jones potential. This is common in particle-based and molecular dynamics simulations that consider infinite potentials at zero distances. As the concentration increases, the probability of finding particles near the wall increases, forming a layered structure. At moderate concentrations, ΦHI = 10%, 15%, particles start to form the first layer next to the wall, and inner particles prefer to stay near the center of the cavity where the steric effects with the particles in the first layer are the weakest. At higher concentrations, ΦHI = 20%, particles form a second layer, since there is not enough space at the center to accommodate them. Thus, the layered structure becomes more pronounced, and the layer separation is determined by the particle size. In the figure, we include the number density of a “bead” system where particles are represented as spheres with an excluded volume with radius a. The spherical cavity has a radius of R = 10a. In this system, particles interact only through far-field hydrodynamics and steric repulsions. However, the layered structure is also observed, becoming a characteristic of highly concentrated confined systems.

FIG. 2.

Number density of the particles within a spherical cavity of radius R = 15 as a function of the particle concentration. (Top) Spheres with rS = 3 and (bottom) cylinders with rC = 2.62 and hC = 2rC. Snapshots for ΦHI = 10% are shown for both systems, while the number density of HI “beads” is included in the inset.

FIG. 2.

Number density of the particles within a spherical cavity of radius R = 15 as a function of the particle concentration. (Top) Spheres with rS = 3 and (bottom) cylinders with rC = 2.62 and hC = 2rC. Snapshots for ΦHI = 10% are shown for both systems, while the number density of HI “beads” is included in the inset.

Close modal

Cylindrical particles exhibit an orientational distribution within the cavity. Similar to liquid crystalline systems,53,54 we define an orientational order parameter λ=123cos2θ1, where cosθ=mnmn, m is the vector parallel to the cylinder’s centerline, and n is the vector pointing from the cavity center to the cylinder’s center-of-mass. A random/disordered structure is characterized by λ = 0, whereas for ordered morphologies, λ = 1 when all cylinders are aligned parallel to the radial direction of the spherical cavity (radial phase), and λ = −1/2 when all cylinders are aligned perpendicular to the radial direction of the spherical cavity (concentric phase). Figure 3 shows the orientational order parameter of cylindrical particles within the cavity as a function of particle concentration. A major result is that for all concentrations, the cylinders are in a disordered state at the center of the cavity and oriented concentrically near the wall. For ΦHI = 5%, the disordered morphology spans all locations in the cavity, whereas for ΦHI = 15%, two concentric regions are observed, separated by a layer where the cylinders are oriented following a ⟨θ⟩ = 35° or 145°. Interestingly, the radial–centripetal ordering, λ = 1, is never observed; we believe that this is due to the small aspect ratio of the cylindrical particles.

FIG. 3.

Orientational order parameter λ of cylindrical particles within a spherical cavity of R = 15. The radius of cylinders is rC = 2.62 and the height is hC = 2rC = 5.24.

FIG. 3.

Orientational order parameter λ of cylindrical particles within a spherical cavity of R = 15. The radius of cylinders is rC = 2.62 and the height is hC = 2rC = 5.24.

Close modal

We now proceed to analyze the short- and long-time diffusive behavior in an attempt to delineate the consequences of long- and short-range HI on the dynamics of the particles. The short-time diffusion coefficient and the mean squared displacement (MSD) of the particles are used to quantify these effects. The particles are suspended in a viscous fluid under zero Reynolds number conditions; they interact with other particles and walls through the HI. Recall that under these conditions, momentum transport is infinitely fast.55–58 Previous studies of confined suspensions have shown that there are multiple factors that originate from HI that should affect the diffusion (mobility) of the particles: (i) the reduction of the particle mobility due to confinement—particles in the bulk diffuse faster than in confined geometries, (ii) the space dependent mobility due to the non-slip conditions at the walls—particle diffusion is zero at the walls, and (iii) the decrease of particle mobility as the concentration increases—there is an interplay between lubrication and long-range HI that becomes important as the average inter-particle distance decreases.29,34,37,59 We seek to quantify and determine the consequences of such factors in a cavity enclosure using our model.

We start by measuring the particles’ radial and tangential short-time diffusivities within the cavity as a function of particle concentration. These transport coefficients are calculated from the relation between MSD and time from Stokes–Einstein following a directional decomposition,60 i.e.,

ΔxR2(ri)=2DR(ri)dt,
(7)
ΔxT2(ri)=4DT(ri)dt,
(8)

for t → 0 and where Δx = x(t + dt)x(t), the radial displacement ΔxR = Δx · x/|x|, the tangential displacement ΔxT = Δx − ΔxR, DR(ri) and DT(ri) are the instantaneous radial and tangential short-time diffusivities at a distance ri from the center of the cavity, and dt is an infinitesimal time interval. DR(ri) and DT(ri) are blocked averaged during a typical simulation at each shell and then averaged over independent simulations. Figure 4 summarizes the diffusion coefficients within the cavity for spherical particles as a function of concentration. In the figure, the diffusivities are normalized by the bulk value, which is defined as single particle diffusivity at infinite dilution.

FIG. 4.

Short-time diffusion coefficients for sphere particles (rS = 3) that are confined in a spherical cavity with R = 15: (top) radial diffusivity and (bottom) tangential diffusivity. The coefficients are normalized by the diffusivity of spherical particles in bulk at infinite dilution, kBT/(6πηrS). The orange dashed line represents the averaged “inner” diffusivities for ΦHI = 10%. The filled shadow area around each curve represents their respective statistical error. The diffusion coefficients for r/a < 2.3 at ΦHI = 20% are missing because there are not enough particles appearing in this zone for the diffusivity measurement due to the layered structure.

FIG. 4.

Short-time diffusion coefficients for sphere particles (rS = 3) that are confined in a spherical cavity with R = 15: (top) radial diffusivity and (bottom) tangential diffusivity. The coefficients are normalized by the diffusivity of spherical particles in bulk at infinite dilution, kBT/(6πηrS). The orange dashed line represents the averaged “inner” diffusivities for ΦHI = 10%. The filled shadow area around each curve represents their respective statistical error. The diffusion coefficients for r/a < 2.3 at ΦHI = 20% are missing because there are not enough particles appearing in this zone for the diffusivity measurement due to the layered structure.

Close modal

We find that the confinement hinders the particle diffusion in both directions. The highest value for the particle diffusivity at ΦHI = 5% is around 60% of the bulk value. The lower particle mobility is directly related to the long-range character of the HI and the non-slip conditions at the walls. As the particle concentration increases, lubrication forces begin to dominate and particle diffusion decreases monotonically. In addition, the short-time coefficients are not constant within the cavity, showing sudden decrease as the particles approach the wall. For ΦHI = 10%, we calculated the averaged “inner” coefficients, which are represented by the orange dashed lines in Fig. 4. Importantly, the decrease in particle mobility at the walls has a stronger effect on the radial diffusion coefficient, indicated by (i) the radial particle mobility decreases at r = 7, when compared with the tangential diffusivity that “feels” the presence of the walls at r = 9, and (ii) the radial diffusivity decreases by 0.2 from r(a) = 7 to the wall contact, whereas the tangential diffusivity only decreases by about 0.12. For a spherical particle near the walls, it is easier to diffuse concentrically than radially. Finally, note that for ΦHI = 20%, the coefficients adopt a non monotonic character, which is correlated with the layered structure of particle density in the cavity. The diffusivity data for r(a) < 2.3 are missing for the case of ϕHI = 20% because the short-time diffusion coefficient at a given location is measured on particles that appear at that location; in other words, the diffusion coefficients are only measured in regions where there is a finite concentration of particles [see Fig. 2 (top)]. Similar observations have also been reported by Aponte-Rivera et al. using a Stokesian dynamics (SD) approach.37 

In Fig. 5, we include the short-time diffusion coefficients of spheres with an excluded volume confined to a spherical cavity, with and without long-range HI (no lubrication forces). Our intention is to change the level of the particle description to isolate the HI contributions. Free draining point-particles do not undergo a space dependent diffusion and the concentration perturbs the diffusion rate only weakly. For HI point particles, the mobility decreases monotonically as the particles approach the walls, and the diffusion rate has a stronger dependence on concentration. Consequently, long-range HI and the zero mobility at the walls are responsible for the non-uniform particle diffusion inside the cavity. Lubrication, on the other hand, imposes a directionality on the short-time mobility—radial diffusion is different than concentric diffusion near the walls—and correlates the particle diffusion with the layered density at high volume fractions. The observed behavior is consistent with the experimental findings.61,62

FIG. 5.

Short-time diffusion coefficients for point-particles (“beads,” rB = 1) that are confined in a spherical cavity with R = 10: radial and tangential diffusivities for HI (left) and free-draining (right) particles. The inset shows the long-time diffusion coefficient as a function of the particle concentration.

FIG. 5.

Short-time diffusion coefficients for point-particles (“beads,” rB = 1) that are confined in a spherical cavity with R = 10: radial and tangential diffusivities for HI (left) and free-draining (right) particles. The inset shows the long-time diffusion coefficient as a function of the particle concentration.

Close modal

It is of interest to compare the short-time diffusion coefficient between spherical and cylindrical particles and to validate the effects of the level of confinement. Figure 6 includes the short-time diffusion coefficients for cylindrical particles (left) in a cavity with R = 15 and for spherical particles (right) in a cavity with R = 30. In the figure, ΦHI = 5% and the results for spherical particles in a cavity with R = 15 are included for reference. Similar to the spherical particles, the cylinders exhibit position-dependent diffusion coefficients, and radial diffusion is affected strongly by the presence of the walls when compared to tangential diffusion. Interestingly, the shape of the cylindrical particles has an important effect on the rate of diffusion. Recall that the cylinders and spheres have the same volume and that the diffusion coefficients are normalized by the bulk diffusivity of spherical particles. Consequently, the short-time diffusion of particles of equal volume is decreased when the geometrical symmetry is broken. Finally, decreasing the level of confinement does not change the qualitative behavior of the short-time mobility but, as the confinement decreases, the inner diffusion coefficients approach the bulk value.

FIG. 6.

Short-time diffusion coefficients for (left) cylindrical particles with rC = 2.62 and hC = 2rC = 5.24 confined in a spherical cavity with R = 15 and (right) spherical particles with rS = 3 confined in a spherical cavity with R = 30. The particle concentration is ΦHI = 5%, and the results for spherical particles with rS = 3 confined in a spherical cavity with R = 15 are included for comparative purposes. The filled shadow area on each curve represents their respective statistical error.

FIG. 6.

Short-time diffusion coefficients for (left) cylindrical particles with rC = 2.62 and hC = 2rC = 5.24 confined in a spherical cavity with R = 15 and (right) spherical particles with rS = 3 confined in a spherical cavity with R = 30. The particle concentration is ΦHI = 5%, and the results for spherical particles with rS = 3 confined in a spherical cavity with R = 15 are included for comparative purposes. The filled shadow area on each curve represents their respective statistical error.

Close modal

Finally, we examine the long time diffusive behavior. We use a generalized Stokes–Einstein relation where the MSD is linearly correlated with a mobility coefficient following a power law by

R(t)R(0)2=Mtα,
(9)

where R is the 3N particle coordinate vector, M is the particle mobility coefficient, and α is the power law exponent that characterizes the type of particle transport. For isotropic diffusion, if α = 1, the particles are diffusive and M = 6D, where D is the particle long-time diffusion coefficient. If α ≠ 1, the transport of Brownian particles is said to be in the anomalous diffusion regime: α < 1 is sub-diffusive, while α > 1 is super-diffusive. For confined systems, as t, the walls impose long-time restrictions on the motion along the confining direction. Therefore, for particles confined in a spherical cavity, the MSD should reach a plateau on a time scale corresponding to the particle diffusion time over the cavity size.

Figure 7 shows the MSD of spherical (top) and cylindrical (bottom) particles in a cavity of R = 15 as a function of particle concentration. The time scale in Fig. 7 is in a diffusion time unit, a2ζ/kBT, where ζ = 6πηa. In the figures, we include the results for ΦHI = 15% of the other particle shape for comparison purposes. The MSDs are collected from ten independent simulations with different random seeds; the particles were able to diffuse more than 300 a diffusion times. As expected, the MSD exhibits diffusive behavior as t → 0 and a plateau when MSD ∼ 200 ∼ R2 at τR, where τR is characteristic particle diffusion time over the cavity. At low concentrations, the diffusive behavior, for spheres and cylinders, spans from t = 0 to t = τR (for 5% spheres, τR ∼ 300 in a2ζ/kBT units). As the particle concentration is increased, the diffusion rate decreases; it is correlated with the short-time behavior, as indicated by the monotonic shift of the MSDs in Fig. 7. For spherical particles at ΦHI ≥ 15%, there is a clear sub-diffusive regime over more than two decades that increases as the concentration is increased. Interestingly, the sub-diffusive regime for cylindrical particles starts at ΦHI = 10%. Note that the MSD for cylindrical particles at ΦHI = 15% is almost equal to the MSD for spherical particles at ΦHI = 20%. In the inset, we have included the MSD of spheres with excluded volume but no lubrication to show how the sub-diffusive regime is never observed when only long-range HI is included, thereby suggesting that the sub-diffusive regime is a lubrication effect. It is then natural to attribute this anomalous diffusion to crowding driven by short-range HI. This regime is then characterized by the diffusion-to-crowding transition τDtC time and by the power law exponent α. In this regime, there are two major features: τDtC and α decrease as the particle concentration increases. In particular, τDtC = 20 and α = 0.85 for cylinders at ΦHI = 10%, while τDtC = 5 and α = 0.65 for cylinders at ΦHI = 15%. After the highly concentrated systems enter the crowding regime and given the fact that there are three dimensional restrictions in their motion as t approaches τR, the particles transition from the slow rate sub-diffusive behavior to the plateau. In  Appendix C, we show that our mean square displacements for spheres agree with the isotropic mean square displacement analysis reported in the literature36,37 both qualitatively and quantitatively. Past work37 reported on the anisotropic behavior of the long-time MSD where both sub-diffusive and subsequent super-diffusive regimes appear at intermediate time scales in the radial direction for highly concentrated systems, while the super-diffusive regime is absent in the tangential direction. In view of that past contribution, our work focuses on the combined effects of different factors on the isotropic long-time MSD instead and provides a first step toward understanding the role of the particle shape on the underlying dynamics.

FIG. 7.

Mean square displacement as a function of time for finite-size particles suspended in a spherical particle of size R = 15. (Top) Spheres with a radius rS = 3 and (bottom) cylinders with rC = 2.62 and hC = 2rC. Inset: the evolution of the mean square displacement for point-particles (“beads”) suspended in a spherical cavity of radius R = 10. In each figure, the MSDs for ΦHI = 15% of the other particle shape are included for comparison purposes. The filled shadow area around each curve represents their respective statistical error.

FIG. 7.

Mean square displacement as a function of time for finite-size particles suspended in a spherical particle of size R = 15. (Top) Spheres with a radius rS = 3 and (bottom) cylinders with rC = 2.62 and hC = 2rC. Inset: the evolution of the mean square displacement for point-particles (“beads”) suspended in a spherical cavity of radius R = 10. In each figure, the MSDs for ΦHI = 15% of the other particle shape are included for comparison purposes. The filled shadow area around each curve represents their respective statistical error.

Close modal

We have used an immersed boundary approach to study the structure and dynamics of suspended spherical and cylindrical particles confined in a spherical cavity.

At low concentrations, the particle number density distribution is uniform in the interior of the cavity. As the concentration increases, a layered structure appears. Cylindrical particles exhibit a random orientation at low concentrations and at the center of the cavity for all concentrations. Excluded volume interactions at the wall force the cylinders to orient concentrically. Interestingly, at high concentrations, the layered morphology of the cylinders is correlated with concentric layers that are separated by a layer of cylinders oriented at an average angle of 35° with respect to the radial direction. Cylinders are never found forming radial morphologies within the spherical cavity.

We used multiple hydrodynamical descriptions of the suspended particles to determine the specific influence of hydrodynamic interactions during the dynamic of the particles (diffusion/mobility): (i) free draining point-particles, (ii) long-range hydrodynamic “beads,” and (iii) finite size particles considering lubrication and long-range HI. We found that long-range HI leads to a position-dependent diffusion of the particles inside the cavity; the particles diffuse faster near the center of the cavity and slower near the walls. The HI also decrease the mobility of the suspended particles when compared with their diffusion in the bulk. The increase in particle concentration also results in a decrease in the particles’ diffusion coefficients; this effect is observed for free-draining and HI particles. However, the concentration decrease in the diffusion rate is stronger when HI are considered. Lubrication forces, or short-range HI, influence the dynamics of highly concentrated suspensions; they generate a direction dependent diffusion, where particles diffuse at a lower rate when moving toward the walls than when moving parallel to the walls. The non-slip conditions at the walls, i.e., zero particle mobility, work synergistically with lubrication forces, resulting in an stronger wall dependence of the diffusion coefficients in the radial direction.

Regarding the long–time dynamics, lubrication gives rise to a sub-diffusive regime at high particle concentrations. The sub-diffusive regime, characterized by the diffusion-to-crowding transition time and the mobility power law exponent, becomes more prominent as the concentration increases.

Introducing cylindrical particles has two major consequences: (i) cylindrical particles have lower short-time diffusion coefficients and (ii) the crowding regime is observed at lower concentrations compared with spheres of equal volume. These observations suggest that the shape of bio-molecules, particles, and polymers could determine their mobility and diffusion inside cells and tissues.

J.L. and X.J. contributed equally to this work.

The development of fast computational codes for simulations of nanoparticles interacting through hydrodynamic and polarization interactions was supported by the Department of Energy, Basic Energy Sciences, Division of Materials Research, through the MICCoM center. The calculations on particle segregation and transport presented in this work were supported by the Department of Energy, Basic Energy Sciences, Division of Materials Research, through the AMEWS center.

There are three important validations that are in order to verify the immersed boundary approach that we used in this work. We start by verifying the fulfillment of Stokes’ law by measuring the sedimentation velocity of a spherical particle that is confined between two parallel walls. Analytical values for this velocity are extracted from previous studies.22,32,33,44,63–65 For the specific IB-pFE-GgEm calculation, the particle radius rS = 5 and the wall distance is H = 15. The periodic boundary condition (PBC) is enforced in the unconfined directions, which are set to a length of 200 (≫H) to avoid the influence of the particle periodic images. The surface of the spherical particle is discretized using 119 surface nodes, resulting in a nodal separation between hmin = 1.334 08 and hmax = 1.905 39 and a smoothing parameter ξIB = 1/0.76hmin. For the calculation, we used a GgEm parameter αGgEm = 0.2 and a mesh resolution with a spacing of 1/2α. The resulting mesh was 60 × 60 × 6 HEX20 elements with 324 886 degrees of freedom. The particle is initially located at (0, 0, d) between the two parallel walls, and it moves parallel to the walls under a sedimenting force with a body force density of (0, 0, 1). The particle’s sedimenting velocity (U) is calculated by averaging particle’s velocities over 100 time steps. Figure 8 (top) shows the sedimenting velocity, normalized by Stokes’ velocity U/U0, as a function of normalized location [(dR)/H]. According to the results, the IB-pFE-GgEm provides an excellent agreement with analytical values.

FIG. 8.

Sphere representation with the immersed boundary method: (top) normalized sedimenting velocity of a spherical particle between two parallel walls as a function of the distance between the particle and the nearest wall. Analytical data are taken from Ref. 44. (Center) Diffusion coefficient of a spherical particle with rS = 3 that is confined in a spherical cavity of size R = 15. SD data are taken from Ref. 37. (Bottom) xy−plane diffusion coefficient for a spherical particle with radius rS = 3 in the center of a slit with heights H = 12, 15, and 20. Analytical data are taken from Ref. 66.

FIG. 8.

Sphere representation with the immersed boundary method: (top) normalized sedimenting velocity of a spherical particle between two parallel walls as a function of the distance between the particle and the nearest wall. Analytical data are taken from Ref. 44. (Center) Diffusion coefficient of a spherical particle with rS = 3 that is confined in a spherical cavity of size R = 15. SD data are taken from Ref. 37. (Bottom) xy−plane diffusion coefficient for a spherical particle with radius rS = 3 in the center of a slit with heights H = 12, 15, and 20. Analytical data are taken from Ref. 66.

Close modal

In addition to Stokes’ law, it is important to verify that our combined approach, between the IB-pFE-GgEm, Fixman’s mid-point algorithm, and the Chebyshev polynomial approximation, is satisfying the fluctuation–dissipation theorem. A sphere that is confined in a spherical cavity and between parallel walls offers a scenario to validate the connection between the diffusion and the fluctuating tensors and the proper calculation of the mobility gradients inside the confined geometries. First, we used three different particle discretizations, using 20, 40, and 60 nodes, for the sphere of size rS = 3 that is confined in a spherical cavity of size R = 15. Figure 8 (center) shows the comparison between the short-time diffusion coefficient of the sphere computed through our IB-pFE-GgEm algorithm and SD algorithm previously reported in the literature.37 These results suggest that even with poor surface descriptions, as long as the IB parameter is appropriately chosen, the IB-pFE-GgEm follows the correct fluctuating short-time behavior. Finally, the long-time diffusion for a sphere confined in a slit as a function of the separation between the walls is shown in Fig. 8 (bottom). Analytical and numerical solutions for this coefficient had been well assessed in the literature.68–71 For the IB-pFE-GgEm, the particle is initially located at the mid-plane between two parallel walls. During the particle’s Brownian movement, its motion is restricted to the plane of symmetry. The surface of the spherical particle is discretized using 20 surface nodes (hmin = 2.198 26, hmax = 2.526 08, and ξIB = 1/0.35hmin). The slit mesh resulted in 60 × 60 × 4 HEX20 elements with 228 872 degrees of freedom. To calculate the error bars of the MSD, five independent simulations for each confinement ratio are performed with a constant time step of 0.002.70 The diffusion coefficients from the IB-pFE-GgEm shows an excellent agreement with analytical and numerical results.

In our immersed boundary description of the finite-size particles, each surface node is linked to its neighboring nodes and the center-of-mass through linear springs with a prescribed stiffness constant k. This parameter controls the stiffness (shape) of particles. If the springs are too weak, particles are deformable and special care must be taken to forbid fluid penetration. On the other hand, if the springs are too stiff, the forces acting on surface nodes will be too large, requiring very small time steps to ensure numerical stability. We performed simulations of spherical particles that are confined in a spherical cavity at high concentrations varying the particle stiffness. We measured the particles’ moment of inertia for each particle ν as follows:

Iν(t)=j=1NSrj(t)2,
(B1)

where Iν(t) is the moment of inertia of particle ν at time t and rj(t) is the distance between node j on surface and the particle center-of-mass. Rigid spheres will have an equal and constant moment of inertia. On the other hand, our semi-rigid particles will show a variation in Iν(t). Figure 9 shows the averaged standard deviation of moment of inertia as a function of time for particles with different k for ΦHI = 20%. The standard deviation is calculated by

σMI(t)=ν=1NIν(t)Iν(0)2N,
(B2)

where Iν(0) is the moment of inertia of particle ν at time 0, which corresponds to a perfect sphere. As one can observe from the figure, in the case that k = 10, the moment of inertia jumps from 0.05 to 0.06 at t ≈ 10, indicating that k = 10 is not stiff enough to maintain particle shape. As the stiffness is increased, the shape variations decrease. In this paper, we used a k = 200 for all simulations.

FIG. 9.

Time evolution of the standard deviation of the sphere moment of inertia as a function of the spring stiffness. The particle volume fraction is ΦHI = 20%, the particle radius is rS = 3, and the number of surface nodes on each particle is 20.

FIG. 9.

Time evolution of the standard deviation of the sphere moment of inertia as a function of the spring stiffness. The particle volume fraction is ΦHI = 20%, the particle radius is rS = 3, and the number of surface nodes on each particle is 20.

Close modal

The equilibrium structure and diffusion in concentrated hydrodynamically interacting spherical particles confined in a spherical cavity were studied using SD by Aponte-Rivera et al.37 There are differences between methods in our and their work such as particles models and HI models. However, one can draw both quantitative and qualitative agreements between results from two methods. Here, we present the comparison of isotropic MSD at different volume fractions [Fig. 15(a) in Ref. 37 and the top panel in Fig. 7 in this work]. Qualitatively, a short-time diffusive regime and a long-time plateau for systems at all concentrations are observed in both works. At intermediate times, a sub-diffusive region emerges for highly concentrated systems only. To make a semi-quantitative comparison, we convert the MSD data in Ref. 37 to the units used in our work at the same level of excluded volume crowding. Specifically, the length and time reported by Aponte-Rivera et al. are normalized by the radius of a sphere with radius rs and the characteristic diffusion time t0 of a sphere with radius rs, respectively (t0=rs2/D0, where D0 = kBT/6πηrs). By multiplying the time, MSD, and volume fraction in the their work by rs3, rs2, and (rs/rs+a)3 (rs = 3a), respectively, we can compare the isotropic MSD directly. Note that the hydrodynamic crowding, ΦHI, and the excluded volume crowding, ΦEV, are equivalent in the work of Aponte-Rivera et al. In contrast, in our work, ΦEV/ΦHI=(rs+a)3/rs3. Thus, the comparison is conducted at the same level of excluded volume crowding but not the same hydrodynamic crowding. For example, the highest concentration in Ref. 37 is 40%, which is comparable to the highest excluded volume crowding in our work, ΦEV = 47% (ΦHI = 20%). For this highest concentrated system, in the work of Aponte-Rivera et al. (after conversion to our units, the converted value is reported in parentheses in the following), the short-time diffusive regime lasts from t ≈ 10−3 (10−2) to ≈10−1 (100), while the MSD increases from 10−3 (10−2) to 10−1 (100); the intermediate sub-diffusive regime lasts from t ≈ 10−1 (100) to t ≈ 101 (102), while the MSD increases from 10−1 (100) to 100 (101); and the MSD reaches a plateau at 101 (102) at t ≈ 102 (103). Good agreement is found between our MSD data and those reported in Ref. 37.

1.
R. J.
Ellis
,
Trends Biochem. Sci.
26
,
597
(
2001
).
2.
R. J.
Ellis
,
Curr. Opin. Struct. Biol.
11
,
114
(
2001
).
3.
H.-X.
Zhou
,
G.
Rivas
, and
A. P.
Minton
,
Annu. Rev. Biophys.
37
,
375
(
2008
).
4.
H.-X.
Zhou
,
Arch. Biochem. Biophys.
469
,
76
(
2008
).
5.
R. J.
Ellis
and
A. P.
Minton
,
Biol. Chem.
387
,
485
(
2006
).
6.
A. P.
Minton
,
J. Pharm. Sci.
94
,
1668
(
2005
).
7.
J. A.
Dix
and
A. S.
Verkman
,
Annu. Rev. Biophys.
37
,
247
(
2008
).
8.
R.
Swaminathan
,
C. P.
Hoang
, and
A. S.
Verkman
,
Biophys. J.
72
,
1900
(
1997
).
9.
B. R.
Terry
,
E. K.
Matthews
, and
J.
Haseloff
,
Biochem. Biophys. Res. Commun.
217
,
21
(
1995
).
10.
M. C.
Konopka
,
I. A.
Shkel
,
S.
Cayley
,
M. T.
Record
, and
J. C.
Weisshaar
,
J. Bacteriol.
188
,
6115
(
2006
).
11.
J.
Skolnick
,
J. Chem. Phys.
145
,
100901
(
2016
).
12.
S. R.
McGuffee
and
A. H.
Elcock
,
J. Am. Chem. Soc.
128
,
12098
(
2006
).
13.
S. R.
McGuffee
and
A. H.
Elcock
,
PLoS Comput. Biol.
6
,
e1000694
(
2010
).
14.
T.
Ando
and
J.
Skolnick
,
Proc. Natl. Acad. Sci. U. S. A.
107
,
18457
(
2010
).
15.
E.
Chow
and
J.
Skolnick
,
Proc. Natl. Acad. Sci. U. S. A.
112
,
14846
(
2015
).
16.
G. B.
Jeffery
,
Proc. R. Soc. London, Ser. A
, Containing Papers of a Mathematical and Physical Character
87
,
109
(
1912
).
17.
M.
Stimson
and
G. B.
Jeffery
,
Proc. R. Soc. London, Ser. A
, Containing Papers of a Mathematical and Physical Character
111
,
110
(
1926
).
18.
J.
Happel
and
R.
Pfeffer
,
AIChE J.
6
,
129
(
1960
).
19.
M.
O’Neill
,
Appl. Sci. Res.
21
,
452
(
1969
).
20.
H. H.
Sherief
,
M. S.
Faltas
, and
S.
El-Sapa
,
Meccanica
52
,
2655
(
2017
).
21.
H.
Lamb
,
Hydrodynamics
(
University Press
,
1924
).
22.
J.
Happel
and
H.
Brenner
,
Low Reynolds Number Hydrodynamics: With Special Applications to Particulate Media
(
Prentice-Hall, Inc.
,
Englewood Cliffs, NJ
,
1965
).
23.
S.
Kim
and
S. J.
Karrila
,
Microhydrodynamics: Principles and Selected Applications
(
Courier Corporation
,
2013
).
24.
G. K.
Batchelor
,
J. Fluid Mech.
44
,
419
(
1970
).
25.
G. B.
Jeffery
,
Proc. London Math. Soc.
s2_14
,
327
(
1915
).
26.
J.
Happel
and
H.
Brenner
,
Low Reynolds Number Hydrodynamics
, (
Springer Science & Business Media
,
1983
).
27.
W. R.
Dean
and
M. E.
O’Neill
,
Mathematika
10
,
13
(
1963
).
28.
M. E.
O’Neill
,
Mathematika
11
,
67
(
1964
).
29.
30.
A. J.
Goldman
,
R. G.
Cox
, and
H.
Brenner
,
Chem. Eng. Sci.
22
,
637
(
1967
).
31.
C.
Pozrikidis
,
J. Fluid Mech.
261
,
199
(
1994
).
32.
P.
Ganatos
,
S.
Weinbaum
, and
R.
Pfeffer
,
J. Fluid Mech.
99
,
739
753
(
1980
).
33.
P.
Ganatos
,
R.
Pfeffer
, and
S.
Weinbaum
,
J. Fluid Mech.
99
,
755
783
(
1980
).
34.
J. W.
Swan
and
J. F.
Brady
,
J. Fluid Mech.
687
,
254
(
2011
).
35.
C.
Oseen
and
N. O.
Zeilon
,
Neuere Methoden und Ergebnisse in der Hydrodynamik
(
Akademische Verlagsgesellschaft m.b.h.
,
1927
).
36.
C.
Aponte-Rivera
and
R. N.
Zia
,
Phys. Rev. Fluids
1
,
023301
(
2016
).
37.
C.
Aponte-Rivera
,
Y.
Su
, and
R. N.
Zia
,
J. Fluid Mech.
836
,
413
450
(
2018
).
38.
X.
Zhao
,
J.
Li
,
X.
Jiang
,
D.
Karpeev
,
O.
Heinonen
,
B.
Smith
,
J. P.
Hernandez-Ortiz
, and
J. J.
de Pablo
,
J. Chem. Phys.
146
,
244114
(
2017
).
39.
C.
Pozrikidis
,
Boundary Integral and Singularity Methods for Linearized Viscous Flow
(
Cambridge University Press
,
Cambridge
,
1992
).
40.
O.
Ladyzhenskaya
,
The Mathematical Theory of Viscous Incompressible Flow
(
Gordon and Beach
,
New York
,
1963
).
41.
H.
Power
and
L. C.
Wrobel
,
Boundary Integral Methods in Fluid Mechanics
(
Computational Mechanics Publications
,
Southampton
,
1995
).
42.
J. F.
Brady
and
G.
Bossis
,
Annu. Rev. Fluid Mech.
20
,
111
(
1988
).
43.
A.
Sierou
and
J. F.
Brady
,
J. Fluid Mech.
448
,
115
(
2001
).
44.
J. W.
Swan
and
J. F.
Brady
,
Phys. Fluids
22
,
103301
(
2010
).
45.
T. A.
Osswald
and
J. P.
Hernández-Ortiz
,
Polymer Processing: Modeling and Simulation
(
Carl Hanser-Verlag
,
Munich
,
2006
).
46.
R.
Cortez
,
SIAM J. Sci. Comput.
23
,
1204
(
2001
).
47.
A.
Kumar
and
M. D.
Graham
,
J. Comput. Phys.
231
,
6682
(
2012
).
48.
49.
P. J.
Atzberger
,
P. R.
Kramer
, and
C. S.
Peskin
,
J. Comput. Phys.
224
,
1255
(
2007
).
50.
H.
Risken
,
The Fokker-Planck Equation
, 2nd ed. (
Springer-Verlag
,
Berlin, Heidelberg
,
1989
).
51.
H. C.
Öttinger
,
Stochastic Processes in Polymeric Fluids
(
Springer-Verlag
,
Berlin, Heidelberg
,
1996
).
52.
P.
Pranay
,
S. G.
Anekal
,
J. P.
Hernandez-Ortiz
, and
M. D.
Graham
,
Phys. Fluids
22
,
123103
(
2010
).
53.
V. N. T.
Tsvetkov
,
Acta Physicicochim. U.R.S.S.
16
,
132
(
1942
).
54.
J. A.
Martínez-González
,
X.
Li
,
M.
Sadati
,
Y.
Zhou
,
R.
Zhang
,
P. F.
Nealey
, and
J. J.
de Pablo
,
Nat. Commun.
8
,
15854
(
2017
).
55.
J. P.
Hernández-Ortiz
,
H.
Ma
,
J. J. d.
Pablo
, and
M. D.
Graham
,
Phys. Fluids
18
,
123101
(
2006
).
56.
J. P.
Hernández-Ortiz
,
J. J. d.
Pablo
, and
M. D.
Graham
,
J. Chem. Phys.
125
,
164906
(
2006
).
57.
J. P.
Hernandez-Ortiz
,
C. G.
Stoltz
, and
M. D.
Graham
,
Phys. Rev. Lett.
95
,
204501
(
2005
).
58.
K. L.
Kounovsky-Shafer
,
J. P.
Hernández-Ortiz
,
K.
Jo
,
T.
Odijk
,
J. J.
de Pablo
, and
D. C.
Schwartz
,
Macromolecules
46
,
8356
(
2013
).
59.
60.
A. E.
Cervantes-Martínez
,
A.
Ramírez-Saito
,
R.
Armenta-Calderón
,
M. A.
Ojeda-López
, and
J. L.
Arauz-Lara
,
Phys. Rev. E
83
,
030402
(
2011
).
61.
M. D.
Carbajal-Tinoco
,
R.
Lopez-Fernandez
, and
J. L.
Arauz-Lara
,
Phys. Rev. Lett.
99
,
138303
(
2007
).
62.
H. B.
Eral
,
J. M.
Oh
,
D.
Van Den Ende
,
F.
Mugele
, and
M. H. G.
Duits
,
Langmuir
26
,
16722
(
2010
).
63.
M. E.
Staben
,
A. Z.
Zinchenko
, and
R. H.
Davis
,
Phys. Fluids
15
,
1711
(
2003
).
64.
R. B.
Jones
,
J. Chem. Phys.
121
,
483
(
2004
).
65.
S.
Bhattacharya
,
J.
Blawzdziewicz
, and
E.
Wajnryb
,
J. Fluid Mech.
541
,
263
292
(
2005
).
66.
M. D.
Corato
,
J.
Slot
,
M.
Hütter
,
G.
D’Avino
,
P.
Maffettone
, and
M.
Hulsen
,
J. Comput. Phys.
316
,
632
(
2016
).
67.
B.
Lin
,
J.
Yu
, and
S. A.
Rice
,
Phys. Rev. E
62
,
3909
(
2000
).
68.
B.
Lin
,
J.
Yu
, and
S. A.
Rice
,
Colloids Surf., A
174
,
121
(
2000
).
69.
E. R.
Dufresne
,
D.
Altman
, and
D. G.
Grier
,
Europhys. Lett.
53
,
264
(
2001
).
70.
N.
Tarantino
,
J.-Y.
Tinevez
,
E. F.
Crowell
,
B.
Boisson
,
R.
Henriques
,
M.
Mhlanga
,
F.
Agou
,
A.
Israël
, and
E.
Laplantine
,
J. Cell Biol.
204
,
231
(
2014
).