Hydrodynamic Enhancement of the Diffusion Rate in the Region

Article pubs.acs.org/JPCB

Hydrodynamic Enhancement of the Diffusion Rate in the Region between Two Fluctuating Membranes in Close Opposition: A Theoretical and Computational Study Martina Pannuzzo,*,† Antonio Grassi,‡ and Antonio Raudino§ †

Department of Computational Biology, Universität Erlangen-Nürnberg, Staudtstrasse 5, 91058, Erlangen, Germany Dipartimento di Scienze del Farmaco and §Dipartimento di Scienze Chimiche, Università di Catania, Viale A. Coria 6, 95125 Catania, Italy

‡

ABSTRACT: Periodic variation of the distance between two weakly adhering bodies gives rise to a huge tangential motions of the sandwiched solvent layer (squeezing flow). Oscillations either can be induced by an external applied field or can spontaneously arise from the coupling with the solvent heat bath. First we calculated by the Navier−Stokes equation the components of the fluid velocity near two oscillating juxtaposed plates. Then we evaluated the influence of plate oscillations on the transport properties of a trace diffusant dissolved at t = 0 in the outer medium for both deterministic and stochastic excitations. By employing both analytical (Fokker−Planck) and coarse-grained molecular dynamics (MD) simulations, we proved that the entry and migration rates of the diffusant sharply increases with the oscillation amplitudes. Enhancement was related to relevant parameters like oscillation frequency, fluid layer thickness, fluid viscosity, and temperature. An extension to the case of oscillating multistacked lamellae has been also made. Theoretical and MD results suggest a significant enhancement of the diffusant flux even in the worse situation of thermally excited small amplitude fluctuations. Excitation arising from other sources (e.g., microwave or ultrasound irradiation of solid−fluid layered systems) could have a dramatic effect on the transport phenomena. Possible implications to relevant biological problems have been discussed.

1. INTRODUCTION AND MOTIVATION Biological systems are built-up by soft objects interacting together through weak forces. The continuous energy exchange with the fluid viscous environment gives rise to endless thermally excited shape (bending) and distance fluctuations that modulate several properties of the living systems. On the other hand, a key event in living systems is the mass transport along the thin aqueous layer confined in the lumen between juxtaposed membrane-bound structures. There are endless examples of biologically relevant layered structures submitted to a heavy traffic of small and large molecules. We mention, for instance, the Golgi apparatus, the Endoplasmic reticulum, the concentric myelin membranes, and the narrow gap between adjacent cells in most tissues. Another interesting topic, occurring on a smaller space and time scales, is the penetration of calcium ions into the narrow gap between a lipid vesicle and a target membrane. Calcium ions efficiently trigger vesicle fusion with a target membrane. The first step of the whole fusion process involves the entry of ions from the external bulk solution into the thin water layer between the tightly adhering membrane-bound bodies. Here they slowly diffuse and tightly bind both to anionic lipids and to fusion proteins (see, e.g., ref 1). A high fusion efficiency requires a fast entry of calcium that, at a first sight, seems to be incompatible with the slowness of hopping diffusive processes2−5 in respect to the higher diffusion coefficient for the free ion.6,7 © 2014 American Chemical Society

Although the two separated issues (namely, membrane oscillations and diffusive transport in confined geometries) have been widely investigated, their mutual coupling and interplay has been generally neglected so far. There are two main different mechanisms that couple distance (or shape) fluctuations with solute transport: (a) a purely hydrodynamic effect related to the oscillationinduced solvent motion near the body boundary (e.g., the membrane−fluid interface) (b) a direct interaction between the oscillating interfaces and the diffusing solute due to long-range intermolecular interactions (e.g., electrostatic forces). In the following we will focus only on the hydrodynamic phenomena described at point (a) neglecting any explicit solute−interface interaction. The diffusion of particles moving inside a fluid within a rapidly changing random environment is an intriguing topic. This complicated effect is usually termed as turbulent or eddy dif f usion. Applications range from oceanographic studies on the mass and temperature transport8,9 stability of flame combustion,10 Received: June 6, 2014 Revised: June 30, 2014 Published: July 3, 2014 8662

dx.doi.org/10.1021/jp505617b | J. Phys. Chem. B 2014, 118, 8662−8672

The Journal of Physical Chemistry B

Article

Figure 1. Cartoon of a disk immersed in an incompressible fluid and oscillating over a flat rigid substrate. Panels A and A* are side (A) and top (A*) views of the system during the compression phase, whereas B and B* are the corresponding drawings during the expansion phase. The black arrows indicate the fluid flow during the expansion/compression motion. Because of fluid incompressibility, the amount of matter entering (leaving) the interplate space is the same during expansion (compression) motion. The red arrows red label the flow of a dilute diffusant moving through an advection-diffusion mechanism, ingoing and outgoing fluxes can be markedly different as described in the text. A real system, closely resembling that one depicted in Figure 1, is shown in Figures 3 and 4.

At a first sight, one could hypothesize that the hydrodynamic fluctuations transport back and forth the diffusant from the bulk fluid to the intermembrane spacing and vice versa. Such a conjecture is valid for the whole mass of the fluid (diffusant + solvent) because of the substantial incompressibility of liquids, but it turns out to be wrong when one is focusing on the fate of the diffusant. Indeed, the periodic squeezing motion affects in a different way the concomitant advection−diffusion process: streamlines entering the narrow gap from the bulk fluid are richer of diffusant molecules, whereas those leaving the gap contain a significantly smaller concentration of diffusants because (a) some of them have been adsorbed onto the reactive lamellae surfaces and (b) in the initial stage the diffusant concentration in the interlamellar gap is lower than in the bulk phase (at t = 0 the concentration in the gap is null), leveling-up to a constant value only at t → ∞. The unbalanced periodic transport of diffusant during the alternating compression/expansion cycles gives rise, after averaging, to an enhanced transport of diffusing particles in the presence of large amplitude fluctuations. The aim of this study is to provide a simplified analytical model for the flux of diffusants that are adsorbed by the reactive cell surface. The models consider both the case of motionless lamellae, where the transport occurs by a pure diffusion mechanism, and that occurring under the hydrodynamic motions induced by lamellae oscillations. In this latter case there is a strong coupling between diffusion and advection that, as we will show later, will enhance the transport rate. Our theoretical predictions will be compared with extensive coarsegrained molecular dynamics (MD) simulations performed on a more realistic molecular model for the interlamellar diffusion. In MD simulations the hydrodynamic effects are considered by explicitly taking into account the chemical nature of the solvent molecules set at close contact with the oscillating lamellae.

dispersion of smoke plumes,11 sediment deposition,12 chromatographic columns,13 and enhancement of the particles mobility in concentrated colloidal suspensions14,15 and in self-replicating living cells,16 just to quote a few exempla. A series of interesting studies has shown over the years that the heat exchange is largely enhanced by turbulent diffusion motions when the heated surfaces vibrate.17−20 This phenomenon explains, at least partially and synergistically with other mechanisms, the ultrasound-induced heat exchange.21 Another potential application of this study is the transdermal drug delivery by oscillating mechanical pressure (sonophoresis22) or small oscillating electric fields (iontophoresis23). Here drugs diffuse faster across the derma when irradiated in the megahertz region. Customarily, the enhancement of drug permeation is ascribed to heating effects or to other ill-defined mechanisms, also in this case a deeper understanding of the basic transport phenomena may help in tailoring more efficient procedures. As we will show later, a likely contribution to the drug transport is related to the convective hydrodynamic motions associated with the cell motions that result strongly enhanced by irradiating the sample by ultrasound or microwaves. In this paper we restrict the analysis to the micron scale by investigating the relative motion of two cells (or tissues madeup of cell layers) surrounded by a thin fluid layer. Cell motion periodically increases and decreases the interlamellar spacing. Because the surrounding fluid is incompressible, what happens is a periodic inner and outer flux of fluid across the thin interstitial fluid film. Next we considered a generic substance (diffusant for short) initially dissolved in the external medium alone. In a later time the diffusant migrates across the thin interlamellar fluid film and is adsorbed (or irreversibly transformed) at the cell surfaces. In the presence of oscillations of the intercellular distance around the equilibrium position, the diffusant motion is strongly affected by hydrodynamic periodic fluxes accompanying the distance oscillations. A cartoon of the considered system is sketched in Figure 1, which gives side (A and B) and top (A* and B*) views of the system during the alternating compression/expansion cycles.

2. THEORY 2.1. Convection−Diffusion Equation in a “Squeezing Flow”. The squeezing phenomenon describes the large flow of a soft material (or a viscous fluid) between two solid surfaces 8663



Article

Once vr has been calculated, we proceed to evaluate vz. By employing the continuity equation together with eq 6, we obtain

approaching each other. During squeezing the gap between surfaces is changing in time and the sample is ejected from the gap. In many applications the gap is small in comparison to the other dimensions of the surfaces, in particular disk diameters (Figure 1), so squeezing flow is mainly associated with thin film hydrodynamics. Its application can be found in various domains like engineering, lubrication, thin films, biofluid dynamics, and rheometry.24 Consider the squeezing motion in a thin film of liquid held between a circular disk of radius R and a parallel flat plate as shown in Figure 1. The distance between the disk and plate is h. We adopt a cylindrical coordinates system (r, z, ψ) with its origin located on the surface of the plate along the symmetry axis of the disk (Figure 1) and impose no-slip boundary conditions on the plate z = 0 vr|z = 0 = vz|z = 0 = 0

∂vz 6U (t ) ⎛⎛ z ⎞ ⎛ z ⎞2 ⎞ 1 ∂ ⎜⎜ ⎟ − ⎜ ⎟ ⎟ =− (rvr ) = ∂z r ∂r h ⎝⎝ h ⎠ ⎝ h ⎠ ⎠

Integrating over z in the interval 0 < z < h subject to vz|z=0 = 0 yields the searched expression for the z-component of the fluid velocity ⎛ ⎛ z ⎞2 ⎛ z ⎞3 ⎞ vz = U (t )⎜3⎜ ⎟ − 2⎜ ⎟ ⎟ ⎝h⎠ ⎠ ⎝ ⎝h⎠

(vr and vz being the r- and z-components of the fluid velocity) whereas on the moving disk located at z = h we have (1b)

vz|z = h = U (t )

(1c)

where U(t) = ∂h/∂t is the velocity of the moving disk along the z axis. In the thin film limit, i.e., h/R ≪ 1, the solution to the Navier−Stokes equation can be obtained via a regular perturbation expansion. To the leading terms, the Navier− Stokes equation reads ∂ 2vr ∂P ∂ ⎛ 1 ∂(rvr ) ⎞ ≈μ ⎜ ⎟+μ 2 ∂r ∂r ⎝ r ∂r ⎠ ∂z

∂ 2vz ∂P ∂ ⎛ 1 ∂(rvz) ⎞ ≈μ ⎜ ⎟+μ 2 ∂z ∂r ⎝ r ∂r ⎠ ∂z

∂C(r,t ) = −∇·(v(r,t ) C(r,t )) + D∇2 C(r,t ) − f (C(r,t )) ∂t

(2a)

(9)

where C(r,t) ≡ C is the local concentration of the diffusant, D is the diffusion coefficient (assumed to be constant), and v(r,t) is the local velocity of the fluid calculated by eqs 6 and 8. The first term in the right-hand side of the transport equation (9) describes the advective term whereas the second (D∇2C) and third ( f(C)) terms account for diffusion and reaction effects, respectively. In the following we adopt cylindrical coordinates r, z, and ψ with 0 < r < R and 0 < z < h (by symmetry, C does not depend on the angular variable ψ). Because we are considering a thin plate configuration where the radii R of the juxtaposed lamellae are much larger than the interlamellar distances h, the space and time variation of the diffusant concentration mainly occurs along the radial coordinate r. Diffusants react at the fluid−plate interfaces where they are adsorbed; however, if the number of binding sites is smaller than that of the diffusant molecules, or if the binding constant is low, the diffusion− advection flux is not significantly changed by the binding of a few diffusant molecules onto the lamellar surface. Therefore, in the following we will set f(C) ≈ 0. Because of the thin layer configuration, we neglect the z-component of the fluid velocity with respect to the r-component vr (they differ by a factor r/h) and vr has been replaced by its averaged value taken over the thickness h of the interlamellar gap

(2b)

where μ is the viscosity coefficient and P is the hydrostatic pressure. For incompressible fluids the continuity equation takes the standard form (1/r)(ϑ/ϑr)rvr + (ϑvz/ϑz) = 0. By symmetry, no angular dependence holds; for that reason, ψdependent terms have been disregarded. The exact solution to eq 2a submitted to the boundary conditions 1a−1c yields vr =

1 ∂P z(z − h) 2μ ∂r

(3)

The radial pressure gradient ∂P/∂r is determined such that continuity is satisfied. Macroscopic continuity requires that the flux out through the wall of the cylinder must balance the flow in through the top

∫0

h

vr ·2πr dz = −πr 2U (t )

(4)

Combining the latter two equations, we find ∂P r = 6μU (t ) 3 ∂r h

(5)

which can be easily integrated in the range 0 < r < R to yield the pressure profile. Combining eqs 3 and 5, eventually we obtain a compact result for the radial component of the fluid velocity 2 r ⎛⎛ z ⎞ ⎛ z ⎞ ⎞ vr = −3U (t ) ⎜⎜ ⎟ − ⎜ ⎟ ⎟ h ⎝⎝ h ⎠ ⎝ h ⎠ ⎠

(8)

Without loss of generality, the oscillation velocity U(t) of the plates can be expressed as U(t) = dh/dt = d/dt (heq + η(t)) where η(t) is the oscillation amplitude of the vibrating plates distance around equilibrium. The instantaneous amplitude η(t) may arise either from deterministic (e.g., by an applied oscillating electric field) or from stochastic forces (e.g., by the coupling of the vibrating plates with the solvent heat bath). This second contribution is ubiquitous and it will be investigated in detail in the next section. 2.2. Solution of the Stochastic Reaction−Advection− Diffusion Equation. Once the fluid velocity components are known (eqs 6 and 8), we proceed to calculate the transport dynamics of a diffusant into the interlamellar gap under the effect of spacing oscillations. The starting point is the advection−diffusion−reaction equation for the diffusant concentration. It takes the usual form

(1a)

vr|z = h = 0

(7)

vr ⇒ ⟨vr ⟩z = − =−

(6) 8664

1 h

1 r ∂η(t ) 2 h ∂t

∫0

h

2 r ⎛⎛ z ⎞ ⎛ z ⎞ ⎞ 3U (t ) ⎜⎜ ⎟ − ⎜ ⎟ ⎟ dz h ⎝⎝ h ⎠ ⎝ h ⎠ ⎠

(10)



Article

By imposing the incompressibility condition given after eq 2, we write down the transport equation (9) into a compact and explicit form ⎛ ∂ 2C ∂η ∂C 1 ∂C ⎞ 1 ∂ = D⎜ 2 + (rC) ⎟+ ∂t r ∂r ⎠ 2h ∂r ∂t ⎝ ∂r

d A m (t ) + [p + G(t )]A m(t ) = 0 dt

where p ≡ and G(t) ≡ (1/2h)(dη(t)/dt). The SDE (15) describes a Langevin equation driven by a multiplicative noise. Such a SDE is common in several fields, ranging from the motional narrowing of the spectral line shape27 to synchronization phenomena28 to noise-induced phase transitions.29 To proceed, we must specify the structure of the random term G(t) related to the stochastic motion of the plates. In the linear approximation, the relative displacement η(t) of two bodies of mass M moving in a medium with a friction coefficient β and connected by a spring of rigidity kel is described by the stochastic equation (d2η(t)/dt2) + γ(dη(t)/dt) + ωo2η(t) = (1/ M)R(t), where ωo2 = kel/M and γ = β/M. The random force R(t) arises from the coupling of the vibrating plates with the heat bath. R(t) is a Gaussian delta-correlated (white) noise with zero mean and autocorrelation function ⟨R(t1) R(t2)⟩ = Γδ(t2 − t1), where Γ is the strength of the noise related to the heat bath temperature through the equipartition theorem. Hence, also η(t) must be a random function with zero mean, ⟨η(t)⟩ = 0, and autocorrelation function

(11)

Equation 11 is a linear stochastic partial differential equation (SPDE) with a time-varying random multiplicative coefficient, ∂η(t)/∂t, to be specified. Its solution must obey proper boundary conditions, specifically: At the center of the of the fluctuating lamellae, r = 0, symmetry imposes

∂C ∂r

=0 (12a)

r=0

At the rim of the fluctuating lamellae, r = R, we assume that diffusant concentration remains constant. That is, those diffusants penetrating the gap are rapidly replaced by others coming from the external solution C|r = R = Cext

(12b)

⎞ kBT ⎛ γ cosh ω t + sinh ω t ⎟e−γt /2 2⎜ * * 2ω Mωo ⎝ ⎠ * 2 kT ≈ B 2 e−ωo t / γ Mωo

Lastly, we assume that at the beginning of the experiment, t = 0, all the diffusants lie outside the interlamellar gap; namely, they are distributed in the external bulk solvent C|t = 0 = 0

⟨η2(t )⟩ =

(12c)

Even if the random term ∂η(t)/∂t is replaced by a deterministic one, the solution to equations like eq 11 is, in general, difficult because the random force is multiplicative and not additive. Often it is convenient to average eq 11 over an oscillation period to get a smooth temporal evolution for the averaged concentration C̅ . This procedure, often called the homogenization or multiple-scale expansion approach, has been developed by various authors over the years.25,26 The basic idea is that of introducing two different times: a long time t defined over the whole time scale of the diffusion process, and a short time τ on the order of plate oscillations. Consequently, the time-varying diffusant concentration are partitioned into two components: C(r,t) = C̅ (r,t) + c(r,τ), where C̅ (r,t) and c(r,τ) are the mean and the fluctuating components of the diffusant concentration, respectively. During the time lapse from t to t + τ, the mean component C̅ (r,t) does not appreciably change, while the time average of the fluctuation part ⟨c(r,t+τ)⟩τ ≡ 1/τ ∫ t+τ t c(r,t+τ) dτ vanishes. However, in the studied case where the r- dependence of the advective term is particularly simple, we were able to find a formally exact analytical solution of eq 11 described by a series of Bessel functions. The proof is reported in the Appendix; here we quote the most relevant formulas C = Cext +

when m ≫ 1

2 1/2

⎞ ⎛ 2 kBT e−ωo / γt ⎟ ⟨A n(t )⟩ = A no exp⎜ −Dλn 2t − 2 2 8Mωo h ⎠ ⎝

(17)

The last step is the calculation of the time-independent coefficients Aom. By exploiting the initial boundary condition 12c, we find (see the Appendix): Aon = −[2Cext/(λnR·J1(λnR))] exp(−kBT/8Mωo2h2). Plugging eq 17 into the averaged eq 13 we find ∞

⟨C⟩/Cext = 1 −

∑ n=1

2J0 (λnr ) λnR ·J1(λnR )

·

⎛ ⎛ ω 2 t ⎞⎞⎞ kBT ⎛ ⎜1 − exp⎜ − o ⎟⎟⎟⎟⎟ exp⎜⎜ −Dλn 2t − 2 2⎜ 8Mωo h ⎝ ⎝ γ ⎠⎠⎠ ⎝

where Jq(λmr) are Bessel functions of order q. Once ⟨C⟩ has been calculated, the local amount Q(r,t) of diffusant molecules that are adsorbed onto the surface of the fluctuating membranes is proportional to the local concentration ⟨C⟩: Q(r,t) ≈ KBind⟨C⟩, where KBind is the binding constant of the diffusant (albeit such an approximation is consistent with the low diffusant concentration assumed throughout the paper, more complex relationships between the fraction of bound molecules and the concentration can be assumed as well. For instance, one could employ the popular Langmuir adsorption equation: Q(r,t) ≈ KBind⟨C⟩/(1 +

(13)

where J0(λmr) is a Bessel function of zeroth order and λm’s are zeros of the transcendental equation J0(λmR) = 0 m− 2.4048 5.5201 , λ2 = , ..., λm = R R R

(16)

with ω* ≡ (1/4γ − ωo ) . This latter approximation holds in the overdamping regime γ/2 ≫ ωo usually found in viscous fluids. Once the statistical properties of the random variable η(t) have been obtained (and, therefore, also those of its time derivative) we easily solve the stochastic eq 15 as sketched in the Appendix; the final result for the time-varying coefficients ⟨An(t)⟩ reads in the overdamped regime 2

(18)

∑ A m(t ) J0 (λmr) m

λ1 =

(15)

Dλm2

1 4

π + O(1/m) (14)

The coefficients Am(t) satisfy the ordinary stochastic differential equation (SDE) 8665



Article

KBind⟨C⟩)). The total amount of bound molecules is found by averaging over the whole surface of the lamellae QTOT(t ) ≈ 4πKBind

∫0

R

⟨C⟩r dr

(19)

a multiplicative factor 2 arises because there are two adsorbing lamellae. In this particular case the integration over r is analytical, whereas if one is using the Langmuir adsorption equation or related relationships, integration must be performed numerically. Exact integration of eq 19 and division by 2πR2 (the area of the two juxtaposed lamellae) yields the wanted expression for the time-varying density ⟨σ(t)⟩ = QTOT(t)/2πR2of the adsorbed diffusant ⎡ ⟨σ(t )⟩ ≈ σmax ⎢1 − ⎢⎣

∞

∑ n=1

Figure 2. Schematic drawing of an array of fluctuating stacked lamellae (side view). Red dots describe a generic diffusant penetrating the interlamellar spacing.

⎛ kBT 4 ·exp⎜ −Dλn 2t − 2 2 λn R ωo 2h2 8 M ⎝

⎞⎤ 2 (1 − e−ωo t / γ )⎟⎥ ⎠⎥⎦

wave vectors q are uniformly distributed over the range qMIN < q < qMAX, where qMIN ≈ π/(N − 1)a and qMAX ≈ π/a, where N is the number of stacked lamellae. Therefore, replacing the sum over q by an integral yields the averaged amount of adsorbed diffusant in the interlamellar space as

(20)

with σMAX ≡ KBindCext and the explicit expression for λn’s is given by eq 14. At small or moderate times eq 20 reduces to ⎡ ⟨σ(t )⟩ ≈ σmax ⎢1 − ⎢⎣

∞

∑ n=1

⎤ 4 2 ⎥ · − λ exp( D t ) eff n ⎥⎦ λn 2R2

⟨⟨σ(t )⟩⟩q = ⟨QTOT(t )⟩q /2πR2 =

(21a)

the effective diffusion constant Deff as modified by hydrodynamic fluctuations reads Deff ≡ D +

(21b)

whereas in the long time limit eq 20 takes a simple form ⟨σ(t )⟩ ∞ ⎤ ⎡ ⎛ kBT ⎞ 4 2 ⎥ ⎢ ⎟ exp( D t ) ≈ σmax 1 − exp⎜ − · − λ ∑ n 2 2 2 2 ⎥⎦ ⎢⎣ ⎝ 8Mωo h ⎠ n = 1 λn R

(21c)

Equations 20−21c are the key results of this paper. It describes the binding kinetics of a diffusant penetrating the spacing between two juxtaposed oscillating lamellae in viscous fluids. 2.3. Extension of the Model to a Fluctuating Multistacks Array. Let us extend the theory developed so far to the rather common situation of a multistacked lamellar array as sketched in Figure 2. Such a situation arises, e.g., in the Golgi apparatus or in the endoplasmic reticulum. At variance with the simple case of two interacting lamellae described above, here the compression−expansion motion of a multistacked array of N lamellae is not characterized by a unique frequency ωo, rather by a distribution of N − 1 different frequencies, ωo → ω(q), the number of allowed frequencies depending on the number of stacks. This property is described by the so-called dispersion relationship that, in the case of identical and equidistant plates, takes the simple form30 ⎛ qa ⎞ ω(q) = ωo · 2 sin⎜ ⎟ ⎝2⎠

∫q

qMAX

MIN

QTOT(t ) dq

(23)

Numerical integration of eq 23 and use of eq 20 enables us to calculate ⟨QTOT(t)⟩q; some results will be reported in section 4. Our result highlights an enhancement of the lateral diffusant transport through lamellar stacks due to broad distribution of the oscillation frequencies in a system of coupled oscillators. 2.4. Numerical Parameters. Different parameters enter in our final equation. A key parameter is the spring constant kel between two membrane-bound bodies in close opposition. Its numerical value depends on the subtle interplay between attractive forces (e.g., ligand−receptor interactions or van der Waals forces) and repulsive interactions (electrostatic, undulation, and hydration forces). Let UTOT(h) be the sum of all attractive and repulsive energy contributions expressed as a function of their mutual distance h. Series expansion around equilibrium distance heq often yields two minima: a deep minimum at very short equilibrium distance (relevant, for instance, in studying fusion events) and a shallow minimum at much longer distances. The spring constant kel is proportional to the compression modulus per unit surface area B for an array of stacked lamellae (an experimental parameter), thus, kel = ∂2UTOT(h)/∂h2|h=heq = πR2B, where πR2 is the surface area of the contacting bodies. Typical values of the compression modulus for lipid bilayers are B ≈ 1012 J·m−4 (ref 31), B ≈ 1014 J·m−4 (ref 32), and B ≈ 1015 J·m−4 (ref 33). This spread in the numerical values of the compression modulus B presents a great variability, dramatically depending on the balance between “stickers” and “repellers”. The temperature was held fixed at 300 K whereas the radius R of the contacting area strongly depends on the geometry of the interacting membranes. In the case of disk-shaped flat interacting structures (e.g., the Golgi apparatus), R is roughly equal to the lateral extension of the interacting bodies. Relating the disk mass M to the radius R, density δ, and disk thickness H, we get M = πR2H·δ. Different results are obtained in the case of weakly interacting spherical cells (or lipid vesicles) because R depends on the cell (vesicle) radius 9 according to34 R3 = const·9 2 .

kBT 8Mγh2λn 2

1 1 2πR2 qMAX − qMIN

(22)

where ωo = kel/M is the maximum frequency (that calculated for two interacting lamellae), q is a wave vector, and a is the interlamellar repeat distance calculated from the lamella’s center of mass (a = h + H, with h the thickness of the fluid layer and H the lamellar thickness, Figure 2). We assume that the 2

8666



Article

Figure 3. (A) Initial setup used in simulations of POPG membranes (cyan) and calcium ions (red). Water molecules are not shown in the snapshot. (B) Frozen POPG membranes (cyan) in the presence of calcium ions (red) mostly adsorbed after 1 μs of simulation. Water molecules are not shown in the snapshot.

small contact area, β is basically twice the friction of two independent spheres calculated by the well-known Stokes relationship: β = 2 × 6πμ9 (μ being the solvent coefficient of viscosity). In pure water, μwater = 10−3 N·s·m−2, but much higher values, strongly depending on the experimental technique and related to the size of the moving particle, have been measured in polymeric solutions or in the cytoplasm.38 It must be borne in mind that in polymeric solvents the viscosity felt by the diffusant and that experienced by a mesoscopic fluctuating body can be markedly different: usually diffusant is much smaller than the polymer gyration radius, so its motion is not significantly affected by the polymer chains.39 On the contrary, the fluctuating body (e.g., a vesicle or cell) is much larger than the polymeric coil size; then friction effects are expected to be huge and related to chains length and concentration.38,40 Hence, by relating the mass M with the spherical cell radius 9 and density δ we find γ ∝ M−7/9 ≈ 107 s−1. On the contrary, for two interacting flat disk-like bodies oscillating perpendicularly to the plane of the disks (Figure 1), β is about twice that found by the Stokes relationship calculated for two independent disks of radius R:41 β = 2 × 16πμR. Relating the disk mass M to their radius R, we find γ ∝ M−1/2 ≈ 107−109 s−1, a higher friction than that found for two interacting spheres of the same mass. Using the calculated values of γ and ωo, we checked that the overdamping regime (γ /2 ≫ ωo) is attained in viscous polymer-containing solvents and/or for weak interlamellar interactions (repeller-filled membranes), a common situation found in biological systems.

Assuming that cells behave as homogeneous elastic bodies, the multiplicative constant is proportional to the adhesion energy per unit surface W and is inversely related to the cell (vesicle) elastic deformation energy E (Young modulus): const ≈ 9/ 2πW(1 − σ2)/E (with σ the Poisson’s modulus). Setting W = 10−4 J·m−2 (a typical value for lipid membrane adhesion in absence of specific stickers31), E ≈ (1−100) × 103 N·m−2 (refs 35 and 36), and σ ≈ 0.3−0.4, we obtain R ≈ 10−7 m for a cell of radius 9 = 10−6 m; that is, the contact region is small in comparison to the particles radius. Recent accurate single-cell measurements gave the mass M and density δ of an isolated cell; typical vales range from 100 (erythrocytes) to 1000 (leucocytes) × 10−15 kg and δ ≈ 1.1 × 1.03 kg·m−3 (ref 37). Making use of the above estimates, we calculated the oscillation frequency of adhering lamellar-shaped and sphere-shaped bodies, respectively. Setting the thickness of the interacting lamellae H = 10−7 m, we evaluated that the oscillation frequency, ωo = (kel/M)1/2, should lie in the range ωo ≈ 108 s−1, a value that does not depend on the mass M (or lateral extension) of the disks. On the contrary, when we consider two interacting spheroidal bodies, in this latter case ωo ≈ 105−106 s−1, a figure slowly decreasing with M (ωo ∝ M−5/8). These results apply to rigid lamellae undergoing only fluctuations in the relative distance. Real lamellae undergo both distance and bending fluctuations; these latter are indirectly included into an effective compression modulus B. Also the intermembrane distance at equilibrium is sensitive to the strength of the attractive and repulsive forces, typically, heq ≈ (2−4) × 10−9 m. The translational diffusion coefficient D ranges from 10−9 m2·s−1 for small-sized molecules or inorganic ions to 10−11 m2·s−1 for proteins. D is related to the diffusant size and solvent viscosity. Likewise, also the friction coefficient γ = β/M of the fluctuating bodies is related to their size and shape: for two identical spherical particles interacting through a

3. MOLECULAR DYNAMICS (MD) SIMULATIONS 3.1. Coarse-Grained Setup. The coarse-grained Martini parametrization42 was used and the following systems were simulated: 8667



Article

Figure 4. (A) Lateral view of two membranes (cyan) free to move along z in the presence of calcium ions (red) after 400 ns. (B) Lateral view of two membranes (cyan) frozen along z in the presence of calcium ions (red) after 400 ns. (C) Calcium contacts with lipid headgroups over the first 1 μs of simulation: in red, the profile corresponding to frozen membranes along z; in black, the profile corresponding to the membranes free to move along z. In the inset is shown the concentration profile over the first 40 ns: at early times the two curves are practically superimposed.

The ion distribution appeared to be stable already after 1 μs (Figure 3B). 3.2. Results of MD Simulations. A strong indication for the oscillation-enhanced transport near a reactive interface can be derived from our molecular dynamics (MD) simulations. We noticed that the binding kinetics of calcium ions to entering in the gap between two adhering anionic lipid membranes behaves differently whenever the membranes are at rest in their equilibrium arrangement or when they undergo thermally excited oscillations. Specifically, as shown in Figure 4, calcium ion migration is faster for oscillating membranes (panel A) than when the membranes are frozen at equilibrium distance (panel B). Although the oscillation-induced diffusion is rather modest in comparison to the frozen situation (panel C of Figure 4) it is worth mentioning that the lipid bilayers form rather stiff pairs (RMSD = 0.07 and 0.7 in the case of frozen and free membranes, respectively) and that the only difference between black and red curves in panel C is the coupling between the lipid bilayers and the heat bath (290 K). Looser bilayer bunches exhibit wider oscillation amplitudes and hence greater ion diffusivity. We want to stress that the MD setup is conceptually similar to the idealized but solvable system reported in Figure 1. This means that, at least qualitatively, analytical and MD methods support each other.

(a) Two flat bilayers made up of 7168 molecules of negatively charged palmitoyl-oleoyl-phosphatidylglycerol (POPG). The bilayers were inserted in a water box (322 568 water beads) in the presence of calcium ions as counterions (3584) placed randomly except in the interlamellar region. The membrane had free borders along the y and x axes to allow for a flux of ions into the intermembrane space (Figure 3A). Position restraints have been applied to lipids along x and y with a force constant of 1 kJ·mol−1·mm−2, whereas they can freely move along z. Periodic boundaries conditions were applied to the three dimensions of the box. (b) The same system as (a) but with position restraints applied also along the z coordinate with a very stiff force constant of 100 kJ·mol−1·mm−2. This choice enables us to simulate and compare flexible (a) and rigid (b) juxtaposed bilayers. Simulations were performed by version 4.6 of the GROMACS simulation package,43 and the MARTINI force field was used for POPG lipids44 and solvent and Ca2+ ions.45 Each system was minimized by steepest descent algorithm after adding each of the molecular species. The Berendsen weak temperature and pressure coupling algorithms46 were utilized with coupling constants of 0.3 and 3.0 ps, respectively. Lipids, water, and ions were separately coupled to a heat bath. The Lennard-Jones potential was smoothly shifted to zero between 9 and 12 Å. For electrostatics, we used the Coulomb potential, which was smoothly shifted to zero between 0 and 12 Å. The essence of coarse-grained (CG) models is to use short-ranged potentials to be computationally efficient. This allowed us to explore cooperative effects of considerable entity like the diffusion of ions from the extralamellar region between the two hypothetical giant vesicles. The simulation time was 2 μs for each of the systems. The time step was set to 0.02 ps.

4. RESULTS AND DISCUSSION Several qualitative (or semiquantitative) conclusions can be drawn from our model and MD simulations. A list follows. 4.1. Hydrodynamic-Enhanced Diffusion. A typical variation of the concentration of a diffusant molecule moving from an outer space (where its concentration is held constant) through a fluid-filled gap between two apposed oscillating bodies calculated by analytical solution (eq 20) of the diffusion−advection equation is shown in Figure 5A,B. In 8668



Article

4.2. Role of the Interparticles Distance. The hydrodynamic enhanced diffusion sharply increases on decreasing the equilibrium thickness h of the thin fluid film squeezed between the oscillating bodies (roughly as exp(|const|/h2)). 4.3. Viscosity Effect. Recalling that in regular fluids the friction coefficient γ is proportional to the fluid viscosity μ and hence is inversely related to the diffusion coefficient D, we conclude that the hydrodynamic-enhanced effective diffusion coefficient defined by the relationship in eq 21b behaves as Deff ∝ μ−1. This is true for simple fluids. In a recent series of papers, we proved that it should be possible to decouple the viscosity effect acting on the diffusant mobility from those affecting the oscillation rate of adhering bodies. This goal is easily reached by using polymer chains soluble in the water phase. It is known that polymers are excluded from the narrow gap between bodies in close opposition by a combination of steric and dielectric effects.39,47 Because of polymer exclusion from the intermembrane gap, the diffusant mobility does not depend on polymer amount. On the contrary, the external viscosity dramatically increases, especially in the case of long polymer chains.38,40 Solvent viscosity strongly damps the oscillation rates of juxtaposed plates. This outcome is very different from that obtained using small-sized viscosity enhancers like sucrose. In the case of small-sized viscosity enhancers, they are homogeneously distributed and the inner and outer solution, so the fluid viscosity reduces both the diffusant mobility and the plates oscillation rate; the net result is a decrease of the diffusant penetration kinetics. Therefore, by employing polymeric viscosity enhancers that selectively modify only the rate of the lamellar motion, we might selectively investigate hydrodynamic effects in transport phenomena. 4.4. Mass and Shape Effects. We investigated the behavior of the hydrodynamic-enhanced diffusion for two relevant geometries: (i) two flat interacting bodies with a large contact area and (ii) two quasi-spherical bodies interacting through a quite narrow contact area. Using the parameters estimated in section 2.4 and recalling that λn ∝ R−1, we get two useful results: (i) When the size (mass) of the flattened interacting bodies is large, the leading mechanism for the interstitial diffusant transport is the hydrodynamic advection term (it decays with the lamellae mass M as M−1/2), whereas the diffusion contribution plays a less relevant role decreasing as M−1. (ii) Slightly different results are obtained for two spherical interacting particles. Here the diffusion term decays as M−4/9, whereas the convective one is proportional to M−2/9. In all cases, the mass of the oscillating body has always an odd effect on the diffusant penetration rate. This latter result is in agreement with the experimental findings: the penetration kinetics of calcium ions through macroscopic arrays of anionic lamellae occurs in the time scale of weeks.48−50 4.5. Diffusion through a Multistack Array. We have proved (eq 23) that the hydrodynamic effects are different in stacks of N oscillating lamellae-like bodies. Such a situation arises, e.g., in the Golgi apparatus or in the endoplasmic reticulum of living cells that are constituted by arrays of submicron-sized membrane-bound flattened sacs51 (cisternae). The hydrodynamic effects on the penetration rate are reported in Figure 6 and have been calculated for N = 2, 10, 20, and 30. The increase of the penetration rate in multistacks with respect to the dimeric case (the lower curve) is evident, the effect rapidly saturating for large values of N. We believe, however, that the calculated transport across a bunch of lamellae is overestimated. Indeed, though the

Figure 5. (A) Variation of the relative averaged concentration ⟨σ(t)⟩/ σMAX of a diffusant penetrating between two oscillating membranes in close opposition against time (s). ⟨σ(t)⟩ has been scaled to the maximum allowable concentration σMAX. At t = 0 diffusant molecules are homogeneously dissolved in the outer solution alone (Figure 3A). Curves have been calculated in the absence (blue line) and in the presence (red line) of hydrodynamic effects. (B) As in panel A. The five red curves have been calculated for increasing frequencies of the interlamellar oscillation. Low frequencies, ωo = 105 s−1, are on the top and high frequencies, ωo = 106 s−1, are on the bottom. At very high frequencies the red curves merge into the blue one calculated by excluding any hydrodynamic effect.

panel A we report the relative diffusant concentration calculated with (red line) and without (blue line) the hydrodynamic effects. In the absence of the squeezing flow, the transport is purely diffusive and it is basically described by an exponential behavior. By contrast, if one introduces hydrodynamic effects associated with periodic compression/dilation of the interstitial fluid film, the diffusant entry is faster in the early stages of the penetration kinetics, then the rate slows down with a long tail. This behavior is common to all the investigated systems, the only difference lying in the magnitude of the hydrodynamic effects. For instance, in panel B of Figure 5 curves have been calculated for increasing values of the interlamellar frequency of oscillation ωo. Notice the dramatic increment of the rate at low frequencies, suggesting that hydrodynamic effects are particularly relevant in loose aggregates. Coarse-grained molecular dynamics calculations reproduce qualitatively similar results (Figure 4). Albeit exactly the same system, the results confirm the validity of our conjecture. 8669



Article

where J0(λmr) is a Bessel function of zeroth order. Moreover, to satisfy the boundary condition 12b, λm’s must be solutions of the transcendental equation J0(λmR) = 0. The coefficients Ao(t) and Am(t) are still unspecified and they will be determined later through the initial condition 12c. The time evolution of the diffusant concentration is obtained by placing expansion A1 into eq 11, yielding ⎡ ∂A m(t )

∞

∑⎢ m=1

⎣

∂t

⎤ + Dλm 2A m(t )⎥J0 (λmr ) ⎦

1 ∂η(t ) = 2h ∂t

⎡ ∂A m(t )

∞

∑⎢ m=1

∂t

m=1

⎝

∂r

⎤ + Dλm 2A m(t )⎥ · ⎦

1 ∂η(t ) 2h ∂t + λm

⎞ + J0 (λmr )⎟ ⎠

(A2)

∂ ∂λm

∫0

∞

∑ A m (t ) · (∫

0

m=1 R

∫0

R

J0 (λmr ) J0 (λnr )r dr

R

J0 (λmr ) J0 (λnr )r dr

J0 (λmr ) J(λnr )r dr )

(A3)

By exploiting the orthogonality relationships ∫ R0 J0(λnr) J0(λmr)r dr = 0 for n ≠ m and ∫ R0 J02 (λnr)r dr = 1/2R2(J12(λnR) + J02(λnR)) = 1/2R2J12(λnR) for n = m, we obtain d A n (t ) ⎡ G d η(t ) ⎤ + ⎢Dλn 2 − n ⎥A n (t ) = 0 ⎣ dt 2h dt ⎦

(A4)

where Gn ≡ 1 + 2λnR[(∂J1(λnR)/∂(λnR))/(J1(λnR)]. By exploiting the differentiation formula of the Bessel functions54 together with the boundary condition J0(λmR) = 0, we find a compact result Gn ≡ 1 + 2λnR =1+

∂J1(λnR )/∂(λnR ) J1(λnR )

2((λnR ) ·J0 (λnR ) − J1(λnR )) J1(λnR )

= −1 (A5)

Exploiting eq A5, the solution to the stochastic differential eq A4 is obtained by exchanging the variable An(t) = e−pt y(t) with p ≡ Dλn2; hence, dy/dt = (−y/2h)(dη/dt). Thus, a formal solution to eq A4 reads A n(t ) = A n(0) exp( −pt −

1 2h

∫0

t

ε(t ′) dt ′)

(A6)

where in the present problem, the Gaussian random variable reads ε(t) ≡ dη(t)/dt ≡ η(t). Starting from the well-known general result55

■

⟨exp[a

APPENDIX. SOLUTION TO THE STOCHASTIC DIFFUSION-ADVECTION EQUATION We look at a solution to eq 11 satisfying the boundary condition 12a of the kind

∫0

t

⎡1 ε(t ′) dt ′]⟩ = exp⎢ a2 ⎣2

t

∫0 ∫0

t

⎤ ⟨ε(t1) ε(t 2)⟩ dt1 dt 2 ⎥ ⎦ (A7)

(here, a = −1/2h). The above result can be greatly simplified when ε(t) ≡ η̇(t), in this particular case ∫ t0 ∫ t0⟨η̇(t1) η̇(t2)⟩ dt1 dt2 = ⟨η2(t)⟩. Hence, the statistical average of eq A6 can be performed once the statistical properties of the stochastic variable η(t) are known. As said in the main text, η(t) has zero mean, whereas its autocorrelation takes the exact expression

∞

∑ A m(t ) J0 (λmr)

⎣

=

assumption of a constant diffusant concentration at the rim of two stacked lamellae is a good approximation in well-stirred solutions, it turns out to be wrong in the case of multistacks. This is because a great number of close-packed rims behaving as sinks for the diffusant molecules lowers the local diffusant concentration (the diffusion equation in the presence of multiple sinks has been widely studied by us in the past52,53). Albeit still rather qualitative (or, at the best, semiquantitative), our analytical and MD results give some first estimates concerning the hydrodynamic effects on the particles mobility across oscillating interstices among arrays of closepacked bodies. Two points are particularly worth mentioning: (a) Hydrodynamic effects (squeezing flow) are greater in multistacked flat bodies, and (b) in the presence of thermal noise alone, the amplification of the hydrodynamic effects is magnified in loose aggregates (namely at low oscillation frequency of weakly bound bodies). These results suggest an explanation for the drug permeation enhancement induced by an oscillating mechanical pressure (sonophoresis22) or by a modulated electric field (iontophoresis23). Indeed, (i) at variance with two juxtaposed bilayers where the resonance occurs over a restricted window of frequency values, in multistacked arrays the resonance involves a broader distribution of different frequencies, and (ii) in irradiated samples, especially on approaching resonance, the exchanged energy between the oscillating bodies and the applied field is much greater than the thermal energy exchanges considered in this paper. Even in the unfavorable case of thermally excited oscillations hydrodynamic effects play a role, as shown by our theoretical and MD simulation. Therefore, it is conceivable that in irradiated systems the wider oscillation amplitudes may induce a much larger enhancement of the diffusion coefficient. Work along these lines is in progress in our laboratory.

m=1

⎛ ∂J (λmr ) 0

∑ A m (t ) · ⎜r

Multiplying both sides of eq A2 by rJo(λnr), integrating over dr in the range 0 to R, and exploiting the identity r(∂Jo(λmr)/ ∂r) = λm(∂Jo(λmr)/∂λm), we get

Figure 6. Variation of the averaged interlamellar diffusant density ⟨⟨σ(t)⟩⟩q (scaled to σMAX) against time (s) calculated for different numbers of stacked lamellae (N = 2, 10, 20, 30).

C = Cext +

∞

(A1) 8670



Article

⎞ −γt /2 kBT ⎛ γ t t cosh sinh e ω + ω ⎜ * * ⎟⎠ 2ω Mωo 2 ⎝ * 2 kT ≈ B 2 e−ωo t / γ Mωo

(8) Piterbarg, L. I.; Ostrowskii, A. G. Advection and Diffusion in Random Media: Implications for Sea Surface Temperature Anomalies; Kluwer: Dordrecht, The Netherlands, 1997. (9) Pasquero, C. Differential eddy diffusion of biogeochemical tracers. Geophys. Res. Lett. 2005, 32, L17603. (10) Echekki, T., Mastorakos, E., Eds. Turbulent Combustion Modeling; Springer: New York, NY, 2011. (11) Gyr, A., Rys, F. S., Eds. Diffusion and Transport of Pollutants in Atmospheric Mesoscale Flow Fields; Kluwer: Dordrecht, The Netherlands, 1995. (12) Graf, W. H. Hydraulics of Sediment Transport; Water Resources Publ.: Highlands Ranch, CO, 1998. (13) Giddings, J. C. ‘Eddy’ Diffusion in Chromatography. Nature 1959, 184, 357−358. (14) Zahn, K.; Mendez-Alcaraz, J. M.; Maret, G. Hydrodynamic Interactions May Enhance the Self-diffusion of Colloidal Particles. Phys. Rev. Lett. 1997, 79, 175−178. (15) Falck, E.; Lahtinen, J. M.; Vattulainen, I.; Ala-Nissila, T. Influence of hydrodynamics on many-particle diffusion in 2D colloidal suspensions. Eur. Phys. J. E 2004, 13, 267−275. (16) Vlad, M. O.; Cavalli-Sforza, L. L.; Ross, J. Enhanced (hydrodynamic) transport induced by population growth in reaction−diffusion systems with application to population genetics. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 10249−10253. (17) Keanini, R. G. Structure and Particle Transport in Second-Order Stokes Flow. Phys. Rev. E 2000, 61, 6606−6620. (18) Carlsson, F.; Sen, M.; Löfdahl, L. Steady streaming due to vibrating walls. Phys. Fuids 2004, 16, 1822. (19) Yu, J. C.; Li, Z. X.; Zhao, T. S. An analytical study of pulsating laminar heat convection in a circular tube with constant heat flux. Int. J. Heat Mass Transfer 2004, 47, 5297−5301. (20) Benavides, E. M. Heat transfer enhancement by using pulsating flows. J. Appl. Phys. 2009, 105, 094907. (21) Legay, M.; Gondrexon, N.; Le Person, S. L.; Boldo, P.; Bontemps, A. Enhancement of Heat Transfer by Ultrasound: Review and Recent Advances. Int. J. Chem. Eng. 2011, http://dx.doi.org/10. 1155/2011/670108. (22) Modi, S. B.; Pokiya, A. B.; Dhandhalya, M. C. A review on Sonophoresis Mediated Transdermal drug delivery system. Am. J. Pharm. Technol. Res. 2012, 3, 257−277. (23) Manabe, E.; Numajiri, S.; Sugibayashi, K.; Morimoto, Y. Analysis of skin permeation-enhancing mechanism of iontophoresis using hydrodynamic pore theory. J. Controlled Release 2000, 66, 149−58. (24) Engmann, J.; Servais, C.; Burbidge, A. Squeeze flow theory and applications to rheometry: A review. J. Non Newtonian Fluid Mech. 2005, 132, 1−27. (25) P.R. Kramer, P. R.; A.J. Majda, A. J. Simplified models for turbulent diffusion: theory, numerical modelling, and physical phenomena. Phys. Rep. 1999, 314, 237−574. (26) Fannjiang, A.; Papanicolaou, G. Convection-enhanced diffusion for random flows. J. Stat. Phys. 1997, 88, 1033−1077. (27) Kubo, R. A stochastic theory of line shape. In Stochastic Processes in Chemical Physics; Shuler, K. E., Ed.; Wiley: New York, NY, 1969; pp 101−127. (28) Cencini, M.; Tessone, C. J.; Torcini, A. Chaotic synchronizations of spatially extended systems as non-equilibrium phase transitions. Chaos 2008, 18, 037125 and references therein.. (29) See, for example: Munoz, M. A. In Advances in Condensed Matter and Statistical Physics; Korutcheva, E., Cuerno, R., Eds.; Nova Science Publ.: Hauppauge, NY, 2004; pp 37−68. (30) See, for example: Goldstein, H.; Poole, C.; Safko, J. Classical Mechanics; Addison-Wesley: Reading, MA, 2001. (31) Bruinsma, R.; Behrisch, A.; Sackmann, E. Adhesive switching of membranes: experiment and theory. Phys. Rev. E 2000, 61, 4253− 4266. (32) Pan, J.; Tristam-Nagle, S.; Kucerka, N.; Nagle, J. Temperature Dependence of Structure, Bending Rigidity and Bilayer Interations of DOPC Bilayers. Biophys. J. 2008, 94, 117−124.

⟨η2(t )⟩ =

(A8)

with ω* ≡ (1/4γ2 − ωo2)1/2, the latter approximate formula being valid in the overdamping limit γ/2 ≫ ωo. Usage of eq A8 leads to ⎛ 1 ⟨A n(t )⟩ = A no exp⎜ −pt + 2 ⎝ 8h ⎞ dt 2 ⎟ ⎠

t

∫0 ∫0

t

•

•

⟨η (t1) η (t 2)⟩ dt1

⎞ ⎛ 2 kBT e−ωo t / γ ⎟ = A no exp⎜ −Dλn 2t + 2 2 8Mωo h ⎠ ⎝

(A9)

Aon

The last step is the evaluation of the coefficients appearing in eq A9. This goal is easily reached by employing the initial condition 12c, C|t=0 = 0, namely, 0 = Cext + ∑∞ n−1⟨An(0)⟩J0(λnr), where ⟨An(0)⟩ is given by eq A9 by setting t = 0. Following the same procedure developed to derive eqs A2 and A3, the coefficients Aon are easily calculated A no = −

⎛ 2Cext kBT ⎞ ⎟ exp⎜ − 2 2 λnR ·J1(λnR ) ⎝ 8Mωo h ⎠

(A10)

Combining eqs A1, A9, and A10, eventually we get the result reported in eq 18 of the main text.

■

AUTHOR INFORMATION

Corresponding Author

*M. Pannuzzo: e-mail, [email protected]. Author Contributions

The authors contributed equally to this work. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Work financially supported by a PRIN 2010/2011, Cod. prog. 2010L96H3K.

■

REFERENCES

(1) McMahon, H. T.; Kozlov, M. M.; Martens, S. Membrane Curvature in Synaptic Vesicle Fusion and Beyond. Cell 2010, 140, 601−605. (2) Albritton, N. L.; Meyer, T.; Stryer, L. Range of messenger action of calcium ion and inositol 1,4,5-trisphosphate. Science 1992, 258, 1812−1815. (3) Al-Baldawi, N. F.; Abercrombie, R. F. Cell Calcium 1995, 17, 422−430. (4) Gabso, M.; Neher, E.; Spira, M. E. Low mobility of the Ca2+ buffers in axons of cultured Aplysia neurons. Neuron 1997, 18, 473− 481. (5) Nakatani, K.; Chen, C.; Koutalos, Y. Calcium diffusion coefficient in rod photoreceptor outer segments. Biophys. J. 2002, 82, 728−739. (6) Li, Y. H.; Gregory, S. Diffusion of ions in sea water and in deep sea sediments. Geochim. Cosmochim. Acta 1974, 38, 703−714. (7) Ribeiro, A. C. F.; et al. Diffusion coefficients and electrical conductivities for calcium chloride aqueous solutions at 298.15 and 310.15 K. Electrochim. Acta 2008, 54, 192−196. 8671



Article

(55) See, for example: Risken, H. The Fokker-Planck equation; Springer: Berlin, 1996.

(33) Salditt, T.; Vogel, M.; Fenzel, W.; Thermal, W. Fluctuations and Positional Correlations in Oriented Lipid Membranes. Phys. Rev. Lett. 2003, 90, 178101. (34) Johnson, K. L. Contact Mechanics; Cambridge University Press: Cambridge, England, 1985. (35) Kuznetsova, T. G.; Starodubtseva, M. N.; Yegorenkov, N. I.; Chizhik, S. A.; Zhdanov, R. I. Neuron Biomechanics Probed by Atomic Force Microscopy. Micron 2007, 38, 824−833. (36) Trickey, W. R.; Baaijens, F. P.; Laursen, T. A.; Alexopoulos, L. G.; Guilak, F. Determination of the Poisson’s ratio of the cell: recovery properties of chondrocytes after release from complete micropipette aspiration. J. Biomech. 2006, 39, 78−87. (37) Grover, W. H.; et al. Measuring single-cell density. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 10992−10996. (38) See, for example: Wirtz, D. Particle-Tracking Microrheology of Living Cells: Principles and Applications. Annu. Rev. Biophys. 2009, 38, 301−326. (39) See, for example: Raudino, A.; Marrink, S. J.; Pannuzzo, M. Anomalous viscosity effect in the early stages of the ion-assisted adhesion/fusion event between lipid bilayers: A theoretical and computational study. J. Chem. Phys. 2013, 138, 234901. (40) See, for example: Wang, P.; Fang, J.; Qin, S.; Kang, Y.; Zhu, D. M. Molecular Weight Dependence of Viscosity and Shear Modulus of Polyethylene Glycol (PEG) Solution Boundary Layers. J. Phys. Chem. C 2009, 113, 13793−13800. (41) Srivastava, D. K.; Srivastava, N. Stokes flow past axially symmetric isolated bodies: a review article. Int. J. Appl. Math. Mech. 2013, 9, 1−53. (42) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The MARTINI force field: coarse grained model for biomolecular simulations. J. Phys. Chem. B 2007, 111, 7812−7824. (43) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (44) Baoukina, S.; Monticelli, L.; Amrein, M.; Tieleman, D. P. The molecular mechanism of monolayer-bilayer transformations of lung surfactant from molecular dynamics simulations. Biophys. J. 2007, 3, 3775−3782. (45) Lee, H.; de Vries, A. H.; Marrink, S. J.; Pastor, R. W. A coarsegrained model for polyethylene oxide and polyethylene glycol: conformation and hydrodynamics. J. Phys. Chem. B 2009, 113, 13186−13194. (46) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684−3690. (47) Raudino, A.; Pannuzzo, M.; Karttunen, M. Combined depletion and electrostatic forces in polymer-induced membrane adhesion: a theoretical model. J. Chem, Phys. 2012, 136, 055101. (48) Feigenson, G. W. On the nature of calcium ion binding between phosphatidylserine lamellae. Biochemistry 1986, 25, 5819−5825. (49) Feigenson, G. W. Calcium ion binding between lipid bilayers: the four-component system of phosphatidylserine, phosphatidylcholine, calcium chloride, and water. Biochemistry 1989, 28, 1270−1278. (50) Swanson, J. E.; Feigenson, G. W. Thermodynamics of mixing of phosphatidylserine phosphatidylcholine from measurements of high affinity calcium binding. Biochemistry 1990, 29, 8291−8297. (51) See, for example: Ladinsky, M. S.; Mastronarde, D. N.; McIntosh, J. R.; Howell, K. E.; Staehelin, M. A. Golgi structure in three dimensions: functional insights from the normal rat kidney cell. J. Cell Biol. 1999, 144, 1135−1149. (52) Baldo, M.; Grassi, A.; Raudino, A. Diffusion-controlled reactions among ligand and receptors clusters. Effects of competition for the ligands. J. Phys. Chem. 1991, 95, 6734−6740. (53) Grassi, A.; Raudino, A. Modeling diffusion-controlled permeation through a fractured surface. Possibility of estimating the averaged distance among the fractures by fluorescence quenching measurements. J. Phys., Chem. 1996, 100, 16934−16939. (54) Gradshteyn, I. S.; Ryzhik, I. M. Tables of Integrals, Series and Products; Elsevier: Burlington, MA, 2007. 8672


Hydrodynamic Enhancement of the Diffusion Rate in the Region

Recommend Documents