Effects of Hydrodynamic Interactions on the Near-Surface Diffusion of Spheroidal Molecules.

We investigated diffusion of spheroidal molecules near a planar surface, accounting for spatially dependent translational and rotational mobilities of molecules resulting from their hydrodynamic interactions with the plane. Rigid-body Brownian dynamics simulations of prolate ellipsoids of revolution of an axial ratio in the range of 1.5 to 3.0, suspended in a viscous fluid, with a no-slip flat boundary confining the suspension were employed. Mobility tensor matrices of molecules were evaluated as functions of spheroids' distance and orientation with respect to the plane. Hydrodynamic interactions with the surface lead to substantial changes of spheroids' translational diffusion coefficients both in the direction perpendicular and parallel to the plane when compared with the values characterizing the bulk diffusion. Moreover, the short-time translational diffusion of molecules, measured in the laboratory frame, both in an unbounded fluid and under the confinement, is non-Gaussian, with much larger deviations from Gaussianity observed in the latter case. In an unbounded fluid, distributions of translational displacements of molecules deviate from those expected for a simple Brownian motion as a result of shape anisotropy. In the presence of the plane, spheroids experience an additional anisotropic drag, and consequently, their mobilities depend on their positions and orientations. Therefore, anomalies in the short-time dynamics observed under confinement can be explained in terms of the so-called diffusing-diffusivity mechanism. Our findings have implications for understanding of a wide range of biological and technological processes that involve diffusion of anisotropic molecules near surfaces of natural and model cell membranes, biosensors and nanosensors, and electrodes.

Introduction

In the biological setting, diffusion of molecules is a major determinant of transport phenomena, signaling, and metabolism.[1−7] It is also a key mechanism of transport in numerous sensing applications.[8−11] The problem of diffusion and adsorption of molecules on surfaces of different origins and of different physicochemical properties, for example, model and natural cell membranes,[12,13] metals,[14,15] minerals,[16,17] and synthetic materials,[18] attracts considerable attention due to its ubiquity and fundamental importance in many research areas such as cell biology, physics and biophysics, chemistry, medicine, and biotechnology.[19−22] Diffusion of proteins toward biological interfaces, such as the cytoplasmic leaflet of the cell membrane, and their interactions with receptors located on those interfaces are important steps in a diverse array of cellular processes, from cell signaling to membrane trafficking.[23] Studies on the mechanism of redox enzymes involve interactions of proteins with surfaces of electrodes.[24,25] In surface-based biosensors, diffusion of target molecules to the sensor may govern the kinetics of binding and ultimately the performance of the sensor.[26] Protein fouling of processing surfaces is of importance in food and pharmaceutical industries,[27] whereas protein adsorption on biomedical implants that are in contact with the bloodstream may lead to formation of clots.[28]

Brownian motions of a body suspended in a viscous fluid near a boundary or a surface are slowed down due to hydrodynamic interactions of the diffusing body with the surface. Moreover, the presence of the boundary introduces anisotropy of Brownian motions of otherwise isotropic particles.[19−32] This phenomenon was described theoretically for spherical molecules, whose translational diffusivities depend on the molecule’s position relative to the boundary,[29,30] and verified in various experiments conducted for nanoparticles.[33−35] For a nonspherical molecule, hydrodynamic interactions with the confining surface result in a mobility tensor matrix, which, even for the simplest case of an axisymmetric shape, is a complicated function not only of the position but also of the orientation of the molecule.[36−40]

Here, we employ the rigid-body Brownian dynamics technique[41−43] to investigate how hydrodynamic interactions with a planar boundary affect local translational dynamics of anisotropic molecules. We consider prolate spheroids of axial ratios falling between 1.5 and 3.0. Such a choice was made because, among the menagerie of known globular proteins, there is an abundance of those whose hydrodynamic shapes can effectively be approximated by an ellipsoid, either triaxial or axisymmetric.[44−46] A few examples among many are hen egg white lysozyme,[47,48] bovine serum albumin,[49] the Fab′ domain of the human IgG class proteins,[50] or bovine pancreatic trypsin inhibitor.[51] Typically, the degree of nonsphericity of such proteins is rather moderate. We should note that, while axial ratios of considered model molecules were derived based on hydrodynamic calculations[52] that we performed for a number of globular proteins, these models lack any protein-specific properties and are essentially no different than, for example, hard spheroidal colloids. Hydrodynamic interactions between diffusing spheroids and the plane are accounted for in simulations by introducing position- and orientation-dependent mobility tensor matrices calculated at each step of the simulation.

We evaluated orientation-averaged translational diffusion coefficients of spheroidal molecules in directions parallel and perpendicular to the flat boundary as functions of the molecule–plane separation and molecule’s size. Diffusion of spheroids in the presence of the surface is drastically slowed down in comparison with that occurring in an unbounded fluid, both in the direction perpendicular and parallel to the plane, in the former case more so than in the latter. Moreover, the short-time translational diffusion of spheroids in an unbounded fluid, measured in the laboratory frame, is non-Gaussian, which is an effect stemming from their anisotropy.[53] A similar effect was previously observed in experiments conducted for uniaxial polymethyl methacrylate ellipsoids of axial ratio 8 that were suspended in water and confined to two dimensions.[54] Variation of spheroids’ mobilities with their position and orientations resulting from the presence of the planar surface further amplifies deviations of translational short-time dynamics of molecules from the simple Brownian picture due to Einstein[55] and von Smoluchowski,[56] which can be explained in terms of the diffusing-diffusivity mechanism.[57,58]

Theory

Near-Wall Mobility Tensor of a Prolate Ellipsoid of Revolution

For an arbitrarily shaped rigid body suspended in an incompressible fluid and influenced by external forces (F) and torques (T), which result in the body’s translational (u) and rotational (ω) velocities, the following linear relation holds, provided that the Reynolds number is low[59] and the creeping flow equations are applicableThe symmetric positive-definite matrix[60] that appears in the above equation represents the resistance (or friction) tensor of the body, Ξ. Its four blocks correspond to the body’s translations (TT) and rotations (RR) and their couplings (TR and RT). The inverse of the resistance tensor matrix results in a symmetric positive-definite matrix of a similar block structure that represents the mobility tensor MThe mobility tensor is related to the diffusion tensor D by the fluctuation-dissipation theorem[61]with kB being the Boltzmann constant and T being the temperature. When an arbitrarily shaped body is suspended near a boundary, such as the planar surface considered here, its resistance and mobility tensor matrices become dependent on the distance and orientation of the body relative to the boundary.[36,62,63] Recently, analytical expressions for the dominant correction to the bulk (i.e., in the unbounded fluid) resistance tensor of an axisymmetric body of a particular shape (such as a rod or a prolate ellipsoid of revolution) due to the presence of a nearby no-slip wall have been derived by Lisicki and coauthors,[37] allowing calculations of position- and orientation-dependent resistance and mobility tensors. These expressions are given in eqs S8–S21 (Supporting Information). Further in the text, we present plots of different components of the mobility tensor matrix of an axisymmetric prolate ellipsoid as functions of the elevation above the surface and the inclination angle. Correction is evaluated using coordinate frames depicted in Figure . The laboratory frame of reference consists of three basis vectors, pointing in Cartesian directions x, y, and z (the basis set consists of three unit vectors {ê, ê, ê}), with the z axis normal to the surface (which is located in the z = 0 plane). The body-wall frame consists of three unit vectors, {v̂1, v̂2, v̂3}, with v̂3 pointing along the long axis of the particle, v̂1 = (ê × v̂3)/ | ê × v̂3|, and v̂2 = v̂3 × v̂1. In the body-wall frame, the resistance and mobility tensors depend on the elevation of the body above the plane (z) and the inclination angle θ, with cos θ = ê · v̂3 (see Figure ).

Figure 1

Left: Schematic representation of a prolate ellipsoid (shown in green) near a planar surface (shown in blue) and of coordinate systems employed. Right: Positions of an ellipsoid may be restricted in simulations to a finite domain (shown in gray) with the plane either present or absent. The latter situation is described in the text as bulk diffusion or diffusion in an unbounded fluid.

The resistance tensor matrix in the body-wall frame {v̂1, v̂2, v̂3} is constructed as a sum of the resistance tensor matrix Ξo of the body suspended in an unbounded fluid and the correcting matrix ΔΞ, whose components, which depend on the shape of the body, are calculated for a given z and θ using analytical expressions that can be found in[37]The bulk resistance tensor Ξo is independent of the position and orientation of the body relative to the plane. Its matrix is diagonaland its components are defined as[64]where η is the fluid viscosity, a is the long semi-axis, and b is the short semi-axis of the ellipsoid. The mobility tensor matrix in the body-wall frame can be calculated by inverting Ξ(z, θ). The mobility tensor in the laboratory frame can be then computed using an appropriate change of basis matrix. For a body residing in the laboratory frame in the xz plane (Figure ), the four blocks of the mobility tensor matrix are of forms[37,38]

Rigid-Body Brownian Dynamics

The following Brownian dynamics propagation scheme for a rigid body, in the absence of external forces, in the laboratory coordinate frame was implemented in-house and employed in the current work[41−43]In the above expression, Δt is the time step, and is the vector describing the position and orientation of the bodyx, y, and z are the positions of the center of the ellipsoid, and Ω1, Ω2, and Ω3 are its orientation. is a random displacement vector during a time step Δt due to the Brownian noisesatisfying the following relations is the position- and orientation-dependent mobility tensor matrix, evaluated in the laboratory frame, and ∇ denotes the divergence operator.

A step of a Brownian dynamics simulation involves updating the three coordinates of the diffusing body’s center of rotation, which is a straightforward operation. Rotational moves are performed using the unbiased protocol described previously by Beard and Schlick,[65] in which a finite rotation Ω = (Ω1, Ω2, Ω3) is represented with the operator UwhereIf we denote some body-fixed vector after and before rotation, respectively, by v and vo, then their relation can be expressed as[66]where the unit vectordefines the axis of rotation. This expression is equivalent to applying the matrix U to rotate vo.[65]

Generation of Brownian motions (eq ) of the body in the presence of the plane requires calculation of translational and orientational divergence terms, which result in correct distributions of body’s positions and orientations at equilibrium.[67] The former are evaluated as standard spatial divergences, whereas the latter are computed using the operator, with v being a body-fixed unit vector.[38,68] In both cases, the random finite difference algorithm (RFD)[38,69] is applied, which, for a position- and orientation-dependent mobility tensor matrix M(r, Ω) (or when orientation is described in terms of the body-fixed vector v, M(r, v)), allows the calculation of translational and orientational divergence terms as[38]In the above expressions, Δs is a Gaussian random variable defined by the momentswhere I3×3 is the identity matrix. Averages in eqs and 25 are computed over a finite set of Δs (or a predefined number of random steps). Some of the tests that we performed to validate the Brownian dynamics algorithm employed in the current work are described in the Supporting Information.

Methods

Parametrization of Spheroids

To derive parameters of spheroidal molecules, we considered several proteins of varying shape anisotropy, whose structures are deposited in the Protein Data Bank.[70] Hydrodynamic properties, that is, sizes and shapes, of these proteins were evaluated directly from their three-dimensional atomic structures using the HullRad package.[52] HullRad uses a convex hull model that accounts for hydration to estimate the hydrodynamic volume of a molecule, calculates the shape factor correction based on a prolate ellipsoid of revolution, and provides the user with parameters such as the hydrodynamic radius, axial ratio, translational and rotational diffusion coefficients, and the intrinsic viscosity of the molecule. Considered proteins, together with their hydrodynamic ellipsoids and axial ratio values, are shown in Figure . Based on results of hydrodynamic calculations, we decided to simulate diffusion of seven ellipsoids. The long semi-axis length of each ellipsoid was set to 73.1Å̊ (which is the length obtained for the most elongated protein from the set presented in Figure , that is, the Zika virus envelope protein (PDB ID: 5JHM)[71]). The length of the short semi-axis was varied between ellipsoids to obtain axial ratio values distributed uniformly in the range of 1.5–3.0.

Figure 2

Exemplary proteins and their hydrodynamic prolate ellipsoids of revolution calculated using the HullRad package.[52] (A) Hen egg white lysozyme (PDB ID: 6LYZ;[72] axial ratio: 1.48). (B) Antibody Fab fragment from Jel 103 (PDB ID: 1MRD;[73] axial ratio: 1.56). (C) F(ab′)2 fragment of the human antibody IgG1 b12 (PDB ID: 1HZH;[74] axial ratio: 1.96). (D) Outer surface protein A from Borrelia burgdorferi (PDB ID: 2OL7;[75] axial ratio: 2.33). (E) Chemokine CCL5/RANTES (PDB ID: 2L9H;[76] axial ratio: 2.40). (F) Chlamydial outer protein N from Chlamydia pneumoniae (PDB ID: 4P3Z;[77] axial ratio: 2.61). (G) Secreted chlamydial protein PGP3 (PDB ID: 4JDM;[78] axial ratio: 2.64). (H) Zika virus envelope protein (PDB ID: 5JHM;[71] axial ratio: 2.82). Drawings were done using the UCSF Chimera package.[79] Proteins are not shown to scale.

Evaluation of Translational Diffusion Coefficients

Translational diffusion coefficients of ellipsoids in directions parallel and perpendicular to the surface were calculated as follows. For an ellipsoid of a given axial ratio and for a given value of the ellipsoid–plane distance in the range from a to 11a, an ensemble of 5 × 106 ellipsoids’ orientations in the laboratory frame was generated using uniformly distributed random rotation matrices[80] (so that the molecules are distributed uniformly over a sphere). Next, for each orientation of the molecule in the generated ensemble, a mobility tensor matrix in the laboratory frame was calculated using analytical expressions given in ref (37), and then all matrices were averaged. Orientation-averaged translational diffusion coefficients in directions x, y, and z of the laboratory coordinate frame, for a given elevation of the ellipsoid above the plane, were calculated from diagonal components of the TT block of the average mobility tensor matrix based on the fluctuation-dissipation theorem given in eq . Our model lacks atomic details of both the molecule and the surface, as well as short-range specific and nonspecific interactions, that would affect the dynamics of the molecule very close to the plane. One may also expect that the continuum description of the solvent breaks for small molecule–surface separations where the granularity of the solvent becomes important.[81] Moreover, as shown in ref (37), the analytical correction to the bulk resistance tensor underperforms for small body-wall separations. Therefore, we decided to consider here only ellipsoid–plane separations above a so that rotations of ellipsoids are not hindered. If, however, one would aim at, for example, studying ordering effects in the presence of planar obstacles,[82] all sterically allowed positions and orientations of molecules should be considered.

Brownian Dynamics Simulations

Short-time diffusion of spheroids, either in the bulk or near the planar surface, was evaluated based on BD simulations using five different setups.

First, for an ellipsoid of a given axial ratio, we generated 3.5 × 106 Brownian dynamics trajectories of 5 × 105 steps each, with Δt (eq ) set to 5 × 10–4. At the beginning of each trajectory, the ellipsoid was so oriented that its long axis was parallel to the z axis of the laboratory coordinate frame. Elevation of the ellipsoid (i.e., the z coordinate of the ellipsoid’s center in the laboratory coordinate frame, Figure ) above the xy plane was chosen randomly, with an equal probability, from the interval (a,5a) (Figure ). During the simulation, the z coordinate of the molecule’s center was restricted to this interval. Restriction on ellipsoid’s translations in the direction normal to the xy plane is imposed using an algorithm described previously[38,83] by rejecting BD steps that result in values of z outside the predefined interval. No restrictions were imposed on coordinates x and y and the molecule’s orientation. Positions and orientations of the ellipsoid along the Brownian dynamics trajectory were collected at different times. As a result, for each considered point in time, we obtained an ensemble of 3.5 × 106 of configurations (positions and orientations) in the laboratory frame. Next, starting from each configuration, a single Brownian dynamics step was performed, with Δt again set to 5 × 10–4. Resulting translational displacements of the ellipsoid in directions x, y, and z of the laboratory coordinate frame, Δx, Δy, and Δz, were collected; their distributions were analyzed as functions of the ensemble evolution time.

Second, for an ellipsoid of a given axial ratio and for the z coordinate of its center fixed at a given value in the range from a to 11a, an ensemble of 108 ellipsoid’s orientations in the laboratory frame was generated using uniformly distributed random rotation matrices.[80] Starting from each configuration, a single Brownian dynamics step was performed, and resulting Δx, Δy, and Δz were collected.

Third, for an ellipsoid of a given axial ratio and for a given fixed value of the inclination angle, in the range from 0° to 90° (with the molecule residing in the xz plane, Figure ), an ensemble of 108 ellipsoid’s positions in the laboratory frame was generated, with the z coordinate of the ellipsoid center chosen randomly, with and equal probability, from the interval (a,5a). Starting from each configuration, a single Brownian dynamics step was performed, and resulting Δx, Δy, and Δz were collected.

Fourth, for an ellipsoid of a given axial ratio, an ensemble of 5 × 108 configurations, with the z coordinate of the ellipsoid center chosen randomly, with an equal probability, from the interval (a,5a), and orientations in the laboratory frame generated using uniformly distributed random rotation matrices, was generated. Starting from each configuration, a single Brownian dynamics step was performed, and resulting Δx, Δy, and Δz were collected.

Fifth, for an ellipsoid of a given axial ratio, a single Brownian dynamics step was performed repeatedly 108 times, starting from a given fixed value of the inclination angle in the range from 0° to 90° and a given fixed value of the z coordinate of the ellipsoid center in the range from a to 11a. Resulting Δx, Δy, and Δz were collected.

Distributions of Brownian Displacements

Distributions of translational displacements (Δx, Δy, Δz) collected from Brownian dynamics simulations of spheroids suspended either in an unbounded fluid or near the surface were analyzed in terms of the second (κ2) and fourth (κ4) cumulants (for the Gaussian distribution, all κ’s for i > 2 equal zero), defined respectively aswith q = x, y, z. Additionally, we calculated the non-Gaussian parameter (or excess kurtosis)[54,57]which approaches zero in the case of the Gaussian distribution and assumes values greater than zero for distributions that are heavy-tailed relative to the Gaussian distribution.

Results and Discussion

Dynamics near the Planar Surface

In Figures –5, we show, using three-dimensional plots and heat maps, the dependence of the mobility tensor matrix of an axisymmetric prolate ellipsoid suspended in a fluid confined by a planar boundary on the position and orientation of the ellipsoid. The mobility tensor matrix is evaluated in the laboratory frame as a function of the molecule’s elevation and inclination angle, as described above, with the long axis of the ellipsoid lying in the xz plane (Figure ). A similar depiction of the mobility tensor matrix of an axisymmetric prolate ellipsoid of the axial ratio 8, under the planar confinement, was presented by De Corato et al.[38] Those authors used the more accurate, numerical approach, namely, the 3D finite element method,[38] whereas in the current work, we evaluate mobility tensors based on analytical expressions for the dominant correction to the bulk friction tensors,[37] as described above. It is thus of interest whether we are able to recover the near-surface dynamics of the ellipsoid described previously. Even more so, as in the original work introducing the correction,[37] the authors compare results of their approach with those of precise multipole simulations only for a case of a rigid rod of the aspect ratio 10 for a set of rod–plane distances and a single value of the inclination angle. Plots presented in the current work are based on calculations performed for an ellipsoid of the axial ratio 4.

Figure 3

Components of the M block of the mobility tensor matrix as functions of the ellipsoid’s elevation above the surface and its inclination angle.

Figure 5

Components of the M block of the mobility tensor matrix as functions of the ellipsoid’s elevation above the surface and its azimuthal angle. (A) Mηa2, (B) Mηa2, (C) Mηa2, and (D) Mηa2.

Translations and Rotations

Nonzero components of the TT block of the mobility tensor matrix (eq ) are shown in Figure . We report values of M elements calculated for the inclination angle θ in the range of 0° to 180° and elevation z in the range of 1 to 10 (in units of a, with a being the long semi-axis of the ellipsoid). Diagonal components M, M, and Mare symmetric around θ = 90°, whereas M, representing the translational coupling in directions x and z (see Figure ), is antisymmetric around θ = 90°.

For any given value of the inclination angle, all diagonal components decrease monotonically with the decreasing elevation, with the motion in the direction perpendicular to the wall (M) being the most affected by hydrodynamic interactions. We note that, for configurations in which the ellipsoid nearly touches the plane, that is, z approaching 1 and θ equal to either 0° or 180°, nonzero mobility values are observed, which is clearly a shortcoming of the analytical approximation as one would expect zero mobility in the limit of the ellipsoid touching the plane.[38] The off-diagonal component Mchanges only slightly with the varying ellipsoid–plane distance. Mequals 0 for θ equal to 0°, 90°, and 180°, regardless of elevation. For other inclination angles, a modest increase in absolute values of Mcan be observed close to the plane. Similarly, an increase with the decreasing body–plane distance was observed in the case of the rod described in ref (37), whereas De Corato and coauthors showed essentially no variation in Mof an ellipsoid for elevation values above a.[38] Overall, however, the off-diagonal M elements assume rather small values.

Nonzero components of the RR block of the mobility tensor matrix (eq ) are presented in Figure . We choose to report M components calculated for the inclination angle in the range of 0° to 180° and elevation in the range of 1 to 4. Similarly, as in the case of M, diagonal components of M are symmetric around θ = 90°, whereas Mrepresenting the rotational coupling in directions x and z is antisymmetric around θ = 90°. For any given value of θ, changes observed in all components of M upon varying the ellipsoid–plane distance between 1 and 4 are only marginal. The inclination angle has only a small influence on M, whereas significant variations with θ are observed for M, M, and M.

Figure 4

Components of the M block of the mobility tensor matrix as functions of the ellipsoid’s elevation above the surface and its azimuthal angle.

Dependencies on θ of all components of M (as well as all components of M), evaluated in the laboratory frame of reference defined in Figure , for large molecule–plane separations converge to those observed in an unbounded fluid. Overall, for the range of ellipsoid–plane distances and inclination angle values considered here, a qualitatively similar behavior has been described in ref (38).

Roto-Translation Couplings

In Figure , elements of the RT block of the mobility tensor (eq ) are depicted. Those exhibit rather complicated behavior close to the surface when θ varies in the range of 0° to 180°.

The Mcomponent (Figure A) that couples the rotation around x to the force applied along y (eqs and 14) is always negative and symmetric around θ = 90°. Therefore, when the ellipsoid is dragged toward y (Figure ), it rotates around x clockwise. The largest absolute values of Mare observed close to the plane for θ approaching either 0° or 180°. In the whole range of ellipsoid–plane distances considered, Mapproaches zero for θ approaching 90°. Rotations around y are coupled with forces acting along x via M(Figure B). Close to the plane, Mis positive for θ approaching either 0° or 180° and negative when θ approaches 90°. This mobility component is also symmetric around 90°. Rotations around y, resulting from the force applied along x, are clockwise for θ, resulting in negative values of M, and counterclockwise otherwise. Rotations around y are coupled with forces acting along z via M(Figure C). This component of the mobility tensor is antisymmetric around 90°, with negative values obtained for θ below 90° and positive values for θ above 90°. Maximal absolute values of Mare observed close to the plane for θ near 45° or 135°. When the force is applied to the center of the ellipsoid lying in the xz plane in the z direction, the ellipsoid rotates around y either clockwise (for θ ∈ (0°,90°)) or counterclockwise(for θ ∈ (90°,180°)). The last nonzero component of M, M(Figure D), couples forces acting along y with rotations around z. Mis antisymmetric around 90°, with positive values obtained for θ ∈ (0°,90°) and negative values above 90°. Consequently, the ellipsoid dragged along y rotates around z either counterclockwise (for θ below 90°) or clockwise (for θ above 90°). All components of M vanish far from the plane (in the considered case of the ellipsoid of the axial ratio 4, they are practically zero when the molecule–plane distance is above 10a), which is consistent with the fact that, in the unbounded fluid, forces acting on the center of the ellipsoid do not induce rotations; that is, there is no roto-translation coupling.

Comparison of the above results, with results of the detailed analysis presented in ref (38), allows us to conclude that, overall, the near-plane dynamics of an axisymmetric prolate ellipsoid resulting from the application of the analytical approximation[37] employed in the current work is, in a qualitative sense, no different than that resulting from finite element calculations.

Translational Drift Velocities

Hydrodynamic interactions between the diffusing spheroidal molecule and the plane result in a quite complex position and orientation dependence of the molecule’s mobility tensor matrix (Figures –5). This, in turn, causes the molecule diffusing near the plane to experience a net drift velocity, which is proportional to the divergence of the mobility tensor matrix.[67,84] Therefore, nonvanishing divergence terms are present in the equation describing the Brownian dynamics algorithm (eq ). In Figure , we show translational drift velocities in directions x () and z () (calculated using the random finite difference method,[38,69]eq ), again assuming that the long axis of the ellipsoid lies in the xz plane (Figure ), as functions of the inclination angle and for four values of the ellipsoid–plane distance. The absence of the translational drift in the y direction stems from the choice of the ellipsoid’s orientation in the laboratory coordinate frame.

Figure 6

Dimensionless drift velocities in directions x (top) and z (bottom) as functions of the ellipsoid’s inclination angle and the ellipsoid–surface distance. Results obtained for an ellipsoid of the axial ratio 4.

For a given distance of the molecule from the plane, displays nonmonotonic character. is negative for θ ∈ (0°,90 ° ) and positive for θ ∈ (90°,180°); it equals zero for inclination angles 0°, 90°, and 180°; and its absolute value decreases with increasing ellipsoid–plane distance. Moreover, for any given distance, the integral of the over θ from 0° to 180° equals zero, which means that translational drift in directions parallel to the plane vanishes, as previously described,[38] when averaged over molecule’s orientations.

The z-component of the translational drift velocity is symmetric around 90°, where it reaches a maximum. Close to the plane, is positive for all values of θ in the range from 0° to 180°. However, when greater ellipsoid–plane distances are considered, z > 2a, displays a different behavior. It initially drops below zero, then increases and becomes positive (for θ equal to ∼15°), and then drops below zero again (for θ equal to ∼165°). This means that, very close to the wall, the molecule experiences the translational drift velocity regardless of orientation, whereas for greater distances, for some values of the inclination angle, the z-component of the velocity vanishes. Absolute values of decrease with increasing distance. More importantly, the z-component of the translational velocity does not have an orientation average of zero.

Our findings regarding the translational drift velocity experienced by an axisymmetric prolate ellipsoid suspended near a plane, while less pronounced, are in a qualitative agreement with those described previously.[38]

Diffusion Coefficients

Recently, diffusion of different-sized proteins (fluorescently labeled antibodies, antibody fragments, and antibody complexes of hydrodynamic radii in the range of ∼29 to ∼237 Å) near surfaces of substrate-supported planar phospholipid bilayers has been studied with the total internal reflection with fluorescence correlation spectroscopy.[12,85] The authors aimed to answer the question whether membranes affect local dynamics of proteins and therefore also the kinetics of binding of proteins to their membrane-associated receptors. They analyzed results of experiments employing the analytical expression for the diffusion coefficient of a spherical particle in the direction normal to the wall proposed by Brenner[30] and concluded that membrane surfaces slow down diffusion of proteins in a size-dependent manner and that the dependence of translational diffusion on the hydrodynamic radius of proteins is stronger than that predicted based on the Stokes–Einstein equation for the diffusion in an unbounded fluid. The authors were able to observe hydrodynamic effects at an evanescent depth length that was several times greater than the hydrodynamic radius of the largest protein considered[12,85] and thus concluded that observed hydrodynamic effects are of long range. We note that diffusion in the direction parallel to membrane surfaces was not considered in the aforementioned studies, and the authors regarded motion of proteins as effectively one-dimensional.

Orientation-averaged diffusion coefficients of spheroidal molecules, in directions parallel (x,y) and perpendicular (z) to the planar surface, evaluated as described above, are given in Figure as functions of the spheroid–plane distance.

Figure 7

Orientation-averaged diffusion coefficients of ellipsoids of different axial ratios as functions of the ellipsoid–surface distance for motions in directions parallel (x, y: indistinguishable) and perpendicular (z) to the surface. Left: dimensionless values; right column: values scaled by isotropic coefficients of an ellipsoid with a similar axial ratio, suspended in an unbounded fluid. Error bars are comparable with symbol sizes. Continuous curves are obtained for a sphere with a radius equivalent to the hydrodynamic radius of an ellipsoid with a given axial ratio, based on the expression given by Brenner.[30]

Brownian motions of molecules are significantly slowed down near the plane with motions in the perpendicular direction being more affected than motions in parallel directions. Close to the boundary, parallel diffusion coefficients of ellipsoids are decreased to 40–50% of values of bulk isotropic diffusion coefficients, whereas perpendicular diffusion coefficients are decreased to 20–30% of bulk values. Molecules with a larger hydrodynamic radius diffuse slower (ellipsoids have equal long axis; therefore, their hydrodynamic radius decreases with the increasing axial ratio). Even for molecule–plane distances that are 10 times greater than the long axis of ellipsoids, values of diffusion coefficients are below those observed in an unbounded fluid.

Hydrodynamic effects on diffusion in considered here model systems are of long range (especially when compared with the effective range of electrostatic interactions that, at physiological salt conditions, is only a few Å), in line with results described in refs (12) and (85), due to the leading inverse dependence of translational diffusion coefficients on the distance from the plane.[37]

Brenner considered a sphere moving toward (or away from) a single planar non-slip surface[30] and derived the following expression for the diffusion coefficient of the sphere in the direction normal to the plane, D(z)whereRH is the hydrodynamic radius of the sphere, Do is its bulk diffusion coefficient, and z is the distance of the sphere from the plane. In Figure , we show numerical values resulting from the Brenner formula applied to spherical molecules of similar hydrodynamic radii as considered ellipsoids. Diffusion of spheres is less affected by the presence of the plane.

In an unbounded fluid, the isotropic diffusion coefficient of a Brownian molecule is linearly dependent on the inverse of the molecule’s hydrodynamic radius. As evidenced in Figure , for any given molecule–plane distance, also orientation- and direction-averaged diffusion coefficients of spheroids suspended near a planar boundary show such a dependence. However, unlike in the case of the bulk diffusion, intercept values of linear fits shown in Figure are below zero and only increase toward zero with the increasing spheroid–plane distance. Slopes of fits also vary with the varying distance from the plane and decrease near the boundary. Moreover, linear dependence on the inverse of the hydrodynamic radius is observed for the orientation-, direction-, and distance-averaged diffusion coefficient.

Figure 8

Orientation- and direction-averaged diffusion coefficient of a spheroidal molecule as a function of molecule’s hydrodynamic radius for different values of the elevation above the surface. For a comparison, a dependence of the diffusion coefficient of a spheroid suspended in an unbounded fluid on the hydrodynamic radius and a similar dependence of the diffusion coefficient averaged over the orientation, direction, and spheroid–surface distance (⟨D⟩) are also shown. In all cases, error bars are comparable with symbol sizes. Continuous lines result from linear regressions.

Non-Gaussian Diffusion of Spheroids in the Bulk and near the Planar Surface

Short-time translational diffusion of spheroids is evaluated in terms of cumulants (eqs and 29) and the excess kurtosis (eq ) that characterize distributions of displacements of the molecule, in directions parallel (Δx, Δy) and perpendicular (Δz) to the plane in a predefined time interval, Δt.

We begin by considering displacements Δx, Δy, and Δz measured at different times on an ensemble of molecules whose positions and orientations evolve in time. In Figure , second cumulants of translational displacements distributions, κ2 (eq ), are shown as functions of the time elapsed from the beginning when all molecules in the ensemble are so oriented that their long axes are parallel to the z axis of the laboratory coordinate frame and their positions above the xy plane are chosen randomly with an equal probability from the interval (a, 5a) (as described above in Section ).

Figure 9

Second cumulants of distributions of displacements of ellipsoids in directions x (similar results, obtained for the y direction, are not shown for clarity) and z as a function of the ensemble evolution time. Top: ellipsoids suspended in an unbounded fluid; bottom: ellipsoids suspended near the planar surface. Error bars are comparable with symbol sizes.

For any given axial ratio of the ellipsoid, the time to reach plateau that is observed in the case of bulk diffusion is somewhat shorter than in the near-plane diffusion case, which is a consequence of the fact that rotational dynamics of molecules is slowed down by hydrodynamic interactions with the plane. This is consistent with depictions of the M block components in Figure showing a modification resulting from the presence of the planar boundary when compared with corresponding bulk values. As ⟨Δx⟩, ⟨Δy⟩, and ⟨Δz⟩ are close to zero, regardless of the elapsed time, that is, the translational drift barely manifests itself in short time intervals during which displacements are measured (in the case of bulk diffusion, average displacements are obviously exactly zero as the divergence of the mobility matrix tensor vanishes), κ2’s are a measure of average squared displacements (eq ) or short-time diffusion coefficients (via eqs and 19). In both cases, diffusion in direction z is initially (when the long axis of molecules is oriented parallel to the z axis of the laboratory coordinate system) faster than the diffusion in direction x (or y). In the case of the bulk diffusion, the short-time diffusion coefficient in direction x increases, the short-time diffusion coefficient in direction z decreases with the ensemble evolution time, and finally, both reach the same value. Orientation-averaged short-time diffusion is thus isotropic. In the case of the near-surface diffusion, a similar behavior of short-time diffusion coefficients for motions parallel and perpendicular to the plane is observed. There is however one important distinction. Short-time diffusion coefficients for motions parallel and perpendicular to the surface converge to different values, with diffusion in the z direction being slower than that in direction x (or y). Short-time orientation-averaged diffusion is thus anisotropic. Ratios of plateau values of κ2 time dependencies, for motions in directions x and z, are in the order of an increasing molecule’s axial ratio (Figure ): 1.45, 1.35, 1.30, and 1.30.

Dependencies of fourth cumulants, κ4 (eq ), of translational displacements distributions on the ensemble evolution time are depicted in Figure . In the bulk, Gaussian diffusion is initially observed, and κ4 vanishes for motions in both x (y) and z directions for short ensemble evolution times.

Figure 10

Fourth cumulant of the distribution of displacements of ellipsoids in directions x (similar results obtained for the y direction are not shown for clarity) and z as a function of the ensemble evolution time. Left: ellipsoids suspended in an unbounded fluid; right: ellipsoids suspended near the planar surface.

As the ensemble evolution time increases and the coherence of orientations of molecules dissipates, values of κ4 also increase, signifying that the short-time diffusion becomes non-Gaussian. For any given value of the molecule’s axial ratio, κ4’s obtained for distributions of Δx (Δy) and Δz converge to the same value. Moreover, κ4’s increase with an increasing axial ratio of the spheroid. In the presence of the plane, short-time diffusion is non-Gaussian even at the initial moment when the coherence of molecules’ orientations is absolute and values of κ4 are then nonzero (Figure ). This is true for both parallel (x, y) and perpendicular motions (z). As the decoherence of orientations progresses in time, κ4’s increase; however, unlike in the case of the bulk diffusion, for any given value of the molecule’s axial ratio, κ4’s of Δx (Δy) and Δz distributions converge to different values. Values of κ4 obtained for Δz distributions are larger than those obtained for distributions of displacements in directions parallel to the plane.

In Figure , non-Gaussian parameters (p, excess kurtosis) for distributions of translational displacements are shown. Clearly, Gaussian diffusion is observed only in the bulk and only for short ensemble evolution times when orientations of spheroids are coherent. As the decoherence progresses, larger p values are observed. In all cases depicted in Figure , excess kurtosis increases monotonically with the progress of the orientational decoherence. Overall, larger p values are observed for molecules diffusing near the plane than for molecules in the bulk. An increase of nearly one order in magnitude is observed for Δx (Δy) displacements and nearly two orders of magnitude in the case of Δz displacements.

Figure 11

Values of the non-Gaussian parameter as functions of the ensemble evolution time for distributions of displacements of ellipsoids suspended either in an unbounded fluid (top) or near the planar surface (bottom) in directions x (similar results obtained for the y direction are not shown for clarity) and z.

In the case of the bulk diffusion, the non-Gaussian parameter converges to zero for short times and increases with an increasing molecule’s axial ratio for moderate and long times. Moreover, plots of p obtained for a given axial ratio for Δx (Δy) and Δz distributions are similar (Figure ).

In the case of the near-plane diffusion, the non-Gaussian parameter decreases with increasing molecule’s axial ratio, for motions in directions parallel to the boundary, for all times. For the perpendicular motion, p decreases with the increasing molecule’s axial ratio for short and long ensemble evolution times. For moderate times, p decreases with increasing axial ratio, which is a consequence of the fact that rates of rotational relaxation of spheroids are different; the faster the rate, the more slender the spheroid.

After a sufficiently long ensemble evolution time, all orientations and positions of a spheroidal molecule in the laboratory frame, within the considered diffusion domain (Figure ) become equally probable, as the distribution over different configurations reaches the stationary Boltzmann distribution. Consequently, time dependencies of p presented in Figure reach plateaus. Let us now consider short-time diffusion evaluated over a subset of molecules at a particular fixed elevation above the xy plane and with orientations distributed uniformly over a sphere. Values of non-Gaussian parameters calculated for corresponding distributions of displacements are given in Figure . In the case of bulk diffusion (Figure , top), p is independent on the position of a spheroid relative to the xy plane. Non-Gaussian parameters characterizing distributions of Δx and Δz (and Δy, not shown in the figure) are equal as orientation-averaged short-time bulk diffusion is isotropic. p increases with increasing molecule’s axial ratio in accordance with data presented in Figure (top). In the case of the near-surface diffusion (Figure , center), p is a decreasing function of the molecule–surface distance. Overall, for a given elevation, distributions of displacements in directions parallel (values of excess kurtosis of Δx and Δy distributions are again equal) to the surface are characterized by smaller p values than the distribution of displacements in the perpendicular direction. In both cases, p depends only weakly on the molecule’s axial ratio. Plateau values observed in Figure are smaller than maximal values of excess kurtosis observed in Figure (center) for the spheroid–surface distance approaching a.

Figure 12

Values of the non-Gaussian parameter calculated for distributions of displacements of molecules in directions x and z. Top: Ellipsoids suspended in an unbounded fluid. Displacements obtained for subsets of configurations in which the z coordinate of molecules’ center is fixed at a given value and their orientations are distributed uniformly over a sphere. Center: Ellipsoids suspended near the planar surface. Displacements obtained subsets of configurations in which the z coordinate of molecules’ center is fixed at a given value and their orientations are distributed uniformly over a sphere. Bottom: Ellipsoids suspended near the planar surface. Displacements obtained for subsets of configurations in which the position of molecules is chosen randomly, with an equal probability, from the interval (a,5a) and their inclination angle θ is fixed at a given value.

We may also consider short-time displacements of molecules near the surface, evaluated over a subset of configurations characterized by a uniform distribution of distance from the plane, with the inclination angle fixed at a particular value. An angular dependency of p is clearly visible in Figure (bottom). As expected, values of p that characterize distributions of displacements at short ensemble evolution times (Figure , bottom) converge to values obtained for θ = 0°.

On the other hand, the short-time diffusion of ellipsoids near the planar surface, evaluated for a particular fixed ellipsoid–boundary distance and for a particular fixed value of the inclination angle, is clearly Gaussian, as evidenced in Figure .

Figure 13

Non-Gaussian parameters calculated for distributions of displacements in directions x, y, and z of ellipsoids suspended near the planar surface. z coordinate of molecules’ center is fixed at 1.25a, and their inclination angle is fixed at a given value.

Exemplary distributions of displacements in directions parallel and perpendicular to the xy plane resulting from Brownian dynamics simulations are shown in Figure (together with Gaussian fits). Apart from the case depicted in Figure C, those have fatter tails than Gaussian distributions.

Figure 14

Distributions of displacements in directions x and z of an ellipsoid of axial ratio 3. (A) Ellipsoid suspended near the planar surface. Displacements obtained for the subset of configurations in which the elevation of the molecule is fixed at 1.25a and its orientations are uniformly distributed over a sphere. (B) Ellipsoid suspended near the planar surface. Displacements obtained for the stationary ensemble of configurations in which the z coordinate of the molecule’s center is chosen randomly, with an equal probability, from the interval (a,5a) and its orientations are uniformly distributed over a sphere. (C) Ellipsoid suspended near the planar surface. The molecule’s elevation is fixed at 1.25a, and its inclination angle is fixed at 30°. (D) Ellipsoid in an unbounded fluid. Displacements obtained for the stationary ensemble of configurations in which the z coordinate of the molecule’s center is chosen randomly, with an equal probability, from the interval (a,5a) and its orientations are uniformly distributed over a sphere.

Non-Gaussianity observed in the case of bulk diffusion originates in anisotropy of molecules (bulk diffusion of isotropic molecules is always Gaussian). In Figure (top), for evolution times close to zero, all molecules in the ensemble have the same orientation with ellipsoids’ axes pointing along axes of the laboratory coordinate frame. Therefore, the mobility tensor matrix of each molecule, evaluated in the laboratory frame, is diagonal. Moreover, it does not depend on the molecule’s position. Translational displacements Δx, Δy, and Δz of each molecule adhere to Gaussian distributionswhere q = x, y, z, and the variance depends on diagonal components of the mobility tensor matrix, (eq ). As the evolution time goes by, the coherence of orientations of molecules vanishes. Also, while translational displacements of a given molecule with a given orientation in the laboratory frame are still Gaussian, distributions that characterize the whole ensemble at a given time, resulting effectively from the mixing of differently oriented molecules whose displacements adhere to Gaussian distributions with different variances (as mobilities in directions of the laboratory coordinate frame axes depend on the orientation of the molecule), are not. Non-Gaussianity increases with the evolution time, with the progressing decoherence. Finally, it reaches a plateau for evolution times that are sufficiently long to allow for the orientational relaxation. Bulk short-time diffusion of anisotropic molecules is thus always non-Gaussian, and largest deviations from Gaussianity are observed for the ensemble of molecules with uniformly distributed orientations, as evidenced in Figure (top). Moreover, increasing values of excess kurtosis are observed for an increasing anisotropy of spheroids. Our findings are in line with results described by Prager,[53] who solved analytically the problem of Brownian motions of axisymmetric molecules that are initially distributed uniformly in a plane and oriented randomly, taking into account the unequal rates of diffusion parallel and perpendicular to the molecular axis of symmetry, and rotational diffusion. Prager concluded that, at times resulting in displacements comparable with sizes of molecules, their diffusion is not strictly Gaussian. Han and coauthors,[54] who investigated Brownian motions of isolated ellipsoids in water confined to two dimensions through digital video microscopy, described that the short-time diffusion of particles is non-Gaussian. Those authors also noted that the short-time non-Gaussian effects should be more pronounced in two than in three dimensions. This is the consequence of the fact that, in three dimensions, the ratio of diffusion coefficients in directions parallel and perpendicular to ellipsoid’s axes saturates at 2 when the ellipsoid’s axial ratio goes to infinity, whereas the ratio of parallel and perpendicular diffusion coefficients saturates at a much larger value under 2D conditions.[54] Indeed, even though non-Gaussian parameters that characterize short-time diffusion in the bulk are distinguishably greater than zero (Figures , top, and 12, top), actual deviations of displacement distributions from Gaussians are rather small (Figure D).

Short-time diffusion of spheroids near the planar boundary shows much larger deviations from the simple Brownian motion than that in an unbounded fluid (Figures –12 and 14). As values of excess kurtosis calculated for different ensemble evolution times are all positive (Figure , bottom), all distributions of translational displacements have tails that decay more slowly than Gaussian. Such a behavior was observed experimentally, for instance, in the case of colloid beads diffusing either along phospholipid bilayers or through entangled F-actin networks[86,87] and nanoparticles (polystyrene) suspended in a solution of polymers (poly(ethylene oxide)),[88] and described theoretically in terms of the diffusing-diffusivity model.[57,58] This model predicts that, in environments in which diffusivity of molecules varies in space, time, or both,[57] the distribution of displacements measured on an ensemble of molecules, in which different molecules have different local diffusivities, becomes non-Gaussian. As we have illustrated above (Figures –5) and as it was described before,[37,38] the mobility of an ellipsoid near a planar surface is a complicated function of its position and orientation relative to the plane. In Figure (bottom), we show that, even at time zero when all molecules in the ensemble are oriented in the same way, values of excess kurtosis characterizing distributions of their displacements are greater than zero, unlike in the case of bulk diffusion. This is a consequence of the fact that positions of spheroids in the ensemble are different and thus their mobilities are different due to distance-dependent hydrodynamic interactions with the plane. As evidenced in Figure , both the distance and orientation dependence of hydrodynamic interactions gives rise to non-Gaussian diffusion of molecules evaluated over the ensemble, while displacements of molecules with a given fixed position and orientation are Gaussian (Figure ). As described previously,[57,58,86] non-Gaussian distributions of short-time displacements evaluated over an ensemble of molecules are generated due to distribution of diffusivities, which in our case, are position- and orientation-dependent (in Figure S4 of the Supporting Information, we show that the distribution of short-time displacements evaluated over an ensemble containing two equal populations of ellipsoids characterized by two different local diffusion coefficients in the (z, θ) space is non-Gaussian, while Gaussian displacements are observed for each population separately). Without such dependence, no excess kurtosis is observed (Figure ). Recently, single particle tracking measurements in a bright-field microscope of spherical colloids diffusing near a planar substrate in gravitational and electrostatic potentials were performed to validate the diffusing-diffusivity mechanism.[57] The authors observed non-Gaussian diffusion in the direction normal to the surface of the substrate; however, they concluded that, in their system, non-Gaussian effects are too small to be seen in displacements parallel to the surface.

Conclusions

We investigated the Brownian motion of spheroidal molecules suspended in a fluid bounded by a planar surface, focusing on effects of nonsphericity of molecules’ shapes and their hydrodynamic interactions with the plane. Confinement introduces anisotropy in motions of spheroids and substantially slows down their diffusion both in the direction parallel and perpendicular to the boundary. Mobilities of molecules depend on their position and orientation relative to the surface, which results in a complicated dynamics near the plane. Distributions of short-time translational displacements of spheroids are non-Gaussian, with tails that decay more slowly than Gaussian, as evidenced by the excess kurtosis analysis. This is a clear deviation from the simple Brownian motion. We observe a similar deviation, although to a much lesser extent, in the case of the short-time diffusion of aspherical molecules in an unbounded fluid. An explanation of these phenomena is presented. Our work provides what we believe is an important insight into a fundamental phenomenon of the Brownian motion of anisotropic molecules under geometric confinement. Moreover, our findings are of interest to a wide audience as diffusion of molecules near various surfaces and boundaries plays a pivotal role in many biological and technological processes.