Kinetics of nanopore transport

12 Pages·0.28 MB·

by  Tom Chou

Kinetics of Nanopore Transport
Tom Chou
DAMTP, University of Cambridge, Cambridge CB3 9EW, ENGLAND
(February 1, 2008)
A nonlinear kinetic exclusion model is used to study osmosis and pressure driven flows through
8 nearly single file pores such as antibiotic channels, aquaporins, zeolites and nanotubules. Two
9 possible maxima in the steady state flux as a function of pore-solvent affinity are found. For small
9 drivingforces, thelinear macroscopic osmotic and hydraulicpermeabilities Pos and Lp,are defined
1 intermsofmicroscopickineticparameters. Thedependencesofthefluxonactivationenergies,pore

n lengthandradius,anddrivingforcesareexploredandArrheniustemperaturedependencesderived.
a Reasonable values for the physical parameters used in the analyses yield transport rates consistent
J with experimental measurements. Experimental consequences and interpretations are examined,

1 and a straightforward extension to osmosis through disordered pores is given.
3
NOMENCLATURE
]
t
f a effective molecular diameter α left entrance rate
o
n total particle number density β right exit rate
s T
. n solute number density γ left exit rate
t s
a k T thermal energy δ left entrance rate
B
m
ℓ section length λ particle mean free path
MFP
- r pore entrance radius p(q) internal right(left) hop rate
p
d
v thermal velocity ξ (α,β,γ,δ,p,q)
n T
o Ap effective pore entrance area χ(0L,R) left(right) solvent mole fraction (0.1)
c E effective activation energies χ(L,R) left(right) solute mole fraction
ξ s
[
Ji(t) flux between sectioniandi+1(1/t) σi site i occupation (0,1)

1 J time average flux(1/t) τ molecular transit times(ξ−1)
i ξ
v L=Nℓ pore length ∆µ solvent chemical potential difference
0
1 L hydraulic permeability ∆P hydrostatic pressure difference
p
0
N number of pore sections ∆Π osmotic pressure difference
0
P osmotic permeability Γ acoustic attenuation coefficient
2 os
X dimensionless driving
0
8
9 I. INTRODUCTION control overallosmotic permeabilities.
/ Mechanisms of fluid transport in confined geometries
t
a Transport through microscopic pores is vital to many are also relevant in industrial applications, especially in
m
biological and industrial processes [1–3]. One phe- separation processes. Both natural and manmade mate-
- nomenon involving flow through microscopic channels is rials such as zeolites contain many molecular-sized, es-
d
n osmosis. Ubiquitous inbiologicalfunctions suchascellu- sentially single-file channels which can selectively absorb
o lar volume control, osmosis and reverse osmosis are also fluids. This size specificity can be exploited in separa-
c exploitedinindustrialprocessessuchassolutionconcen- tion of a mixture of linear and branched chain alkanes,
v: tration, filtration, and catalysis. where the zeolite acts like a sponge absorbing only the
i Living cells must transport nutrients, ions and water desired specie(s) [3]. Confining particles in zeolite pores
X across their lipid membranes. Various biological struc- canalsoservetocatalyzereactions[5,6]. Therefore,there
r turesandmachinery,suchasactivetransporters,cotrans- hasbeenmuchresearchintofluidstructureandtransport
a
porters, antiporters, and simple channels or pores have within confined pores [7–16].
evolved to perform these tasks [1,2]. Of central impor- Although the equilibrium properties of osmotic “pres-
tance in living organisms is the flow of water. Certain sure” are understood in terms of macroscopic forces [17]
cells, such as nephrons, have a permeability to osmot- and statistical mechanics [18], nonequilibrium osmotic
ically driven water transport too high to be explained transport has been less well categorized. One approach
by simple permeation through an oily lipid membrane uses macroscopic flow equations (such as the parabolic
bilayer. In addition to generic antibiotic pores such as Poiseuillean fluid velocity profile through a pipe [19]);
gramicidin, pores which are apparently water specific these,appliedtostructuresofmoleculardimensions lead
(aquaporin-CHIP) have been identified [4]. By inserting to inaccurate prediction of membrane pore density, or
aquaporins into their membranes, cells can in principle pore radius, particularly in biological examples [2,20].
1
Use of macroscopic parameters such as viscosity, or hy- ate the system-dependent parameters used in this work,
draulic flow, to microscopic transport has led to mod- provided accurate force fields are known.
els that do not accurately describe most of the relevant In the next Section, we review the usual linear
microscopic physical processes [21,22,20] and many con- phenomenological expressions used to describe osmosis
flicting viewpoints [2]. A recent review of the historical or pressure driven flows through semipermeable mem-
controversies and misconceptions is given by Guell and branes and identify all macroscopic variables. A one-
Brenner [17], who present a clarifying macroscopic de- dimensional symmetric exclusion model is then formu-
scription of osmosis. Finkelstein [2] also mentions the lated for particle dynamics within a pore that connects
many conceptual disagreements associated with osmosis two fluid reservoirs. The rate coefficients used are ap-
and nonequilibrium thermodynamics, particularly when proximated using kinetic theory near thermodynamic
applied to water transport across biomembranes. equilibrium.
In this paper, we consider molecularly-sized pore that In the Results and Discussion, we explore the model
are nearly single-file. Many natural examples of pores as a function of reasonable parameters and present the
are of molecular sizes. Zeolites typically have pores with mean flow rates in a series of plots. The macroscopic
radii∼2−10˚Aandlengthfromafewnanometerstomil- variables typically used are then related to the rate con-
limeters. BiologicalchannelsIntegralmembraneproteins stantsusedinthemicroscopictheory. Nonlinearitiesand
that form channels in biological systems are difficult to Arrenhius temperature dependences are also discussed.
prepare for X-ray crystallography, accurate dimensions In the Conclusions, we summarize our results and intro-
of membrane water channels are not yet available. How- duceextensionstobestudiedinthefollowingpaper. The
ever,spanningcellmembranes,theyare∼50˚Ainlength. Appendix gives a simple result for pores with disordered
Based on sequential analogies with antibiotic pores and interiors, or hopping rates.
electroncryocrystallography[26,27]biopores such as the
water transporting aquaporins are believed to have pore
diameters ∼2−4˚A. Thus, water movement inside these II. MEMBRANE AND PORE TRANSPORT
pores is statistically single-file, with “overtaking” rarely MODELS
if at all occuring.
A. Macroscopic Phenomenological Expressions
(L) (R)
Here, we briefly review the macroscopic descriptions
which hitherto have typically been used to model osmo-
L=N l
sis. The typical starting assumption is that the solvent
R
α q β P flux is a function of the solventchemicalpotentialdiffer-
p
ence between two reservoirs (see Figure 1) that do not
PL γ δ
exchange solute. For two reservoirs with solvent at dif-
ferent equilibrium chemical potentials,
J =J(∆µ)= c (∆µ )k = d (−1)k(∆χ )k,
k 0 k 0
k=1 k=1
X X
FIG. 1. Schematic of osmosis and pressure driven flow (2.1)
through membrane pores. The reservoirs (L) and (R) are
assumed infinite. The rates ξ = (α,β,γ,δ,p,q) are condi- wherethelastequalityassumesTL =TR, PL =PR,and
tional solvent entrance probabilities at the ends of the pore c andd areexpansioncoefficientswhichmaydependon
k k
indicated (see text). For typical experiments, the membrane external parameters such as temperature, pore-particle
is impermeable to solute, ∆P =PR−PL, (L) contains pure interactions, etc. They may also depend on the total
solvent, and J indicates thedirection of solvent flux. soluteconcentrationinthe reservoirs. Equation(2.1)as-
sumes linear response and expands ∆µ in powers of sol-
Using a symmetric exclusion model [28], thermody- vent density (a virial expansion). The microscopic inter-
namics, and simple kinetic theory, we formulate an ap- actions are nearly identical for solventin either reservoir
proximate model for transport through nearly single-file (nearly allvirialterms for the solventchemicalpotential
pores [16] which can be solved exactly. Estimates of the cancel) except for the infrequent encounters with solute
microscopic kinetic parameters used allow the theory to in the solution(R) side. To lowestorder,the (low)solute
be general enough to predict trends and to explore de- concentration difference gives rise to the free energy dif-
pendences on physical parameters such as pore geome- ference which can be seen by using χ = 1−χ , such
s 0
try, microscopic interactions, and bulk fluid properties. that J can be written in terms of solute mole fraction
Our model complements, although is more flexible than difference∆χ ≡χR−χL. Asimilarexpansioninhydro-
s s s
detailed numerical simulations. Microscopic motions de- static pressure gives van’t Hoff’s result that the solvent
rived from simulations can in principle be used to evalu-
2
flowinducedbyanimpermeablesolute isthermodynam- t contains a particle (σ (t)=1) or is empty (σ (t)=0).
i i
ically balanced (J = 0) by a hydrostatic pressure differ- The discrete dynamics are defined by the follow rules:
ence equal to ∆p = −∆Π = −k T(nRχR −nLχL) ≃ Duringanyinfinitesimaltimeintervaldt,randomlychose
B T s T s
−k Tn ∆χ , where n ≃ nL ≃ nR are the total par- a particle with label 0 ≤ i ≤ N +1, where i = 0,N +1
B T s T T T
ticle number densities. Small currents between (R) and corresponds to choosing the (L),(R) reservoir. If i = 0
(L)arethus assumed(fairly accurately)tobe linearJ ≈ (the left reservoir) is chosen during dt, and σ (t) = 0,
1
L (∆Π−∆p). Under isobaric conditions, and χ ≪ 1, thenwefillthefirstsite(σ =0)→(σ (t+dt)=1)with
p s 1 1
osmotic flow is often represented by the osmotic perme- probability αdt. If i = N +1 (the right reservoir) and
ability P (also denoted P ) where J ≃P (nR−nL). σ (t) = 0, then we extract a particle from (R) and set
os f os s s N
Theconstantsofproportionality,L andP determine σ (t+dt) = 1 with probability δdt. These steps allow
p os N
the rates of flow and have been the subject of consider- injection of the end sites of the pore from the reservoirs
able attention. Measurements of membrane permeabili- if and only if the end sites are empty at time t. If a par-
ties and the number of pores per area, can yield single ticleati=1ischosen,weemptyitintothe(L)reservoir
pore conductivities. Conversely, the number of pores in with probability γdt, and if site i=2 is empty, we move
a cell or vesicle can also be determined if single pore the particle to the right with probability pdt. Similarly,
conductances can be accurately measured or modelled. if a particle at site i = N is chosen, it moves into the
Fully microscopic approaches,such as molecular dynam- (R) reservoir with probability βdt, and if σ (t) = 0,
N−1
ics (MD) simulations reveal detailed microscopic infor- it moves to the left with probability qdt. If an interior
mation on how fluids particles adsorb to the pore inte- site, 1 < i < N, is chosen, and i+1(i−1) is empty,
rior [23,24]; however, the time scales readily achievable then we move the particle to the right(left) with prob-
by MD (∼ ns) cannotreadily measuresteady state flows ability pdt(qdt) . The model assumes that p and q are
throughalongpore,orallowexplorationoflargeparam- constantalongalltheinteriorsites: Effectsofsitedepen-
eter dependent flowregimes. To achievechannelflowsin dentp=q are derivedinAppendix A.All the movesare
MDsimulations,anexternalpotential(suchasgravityor allowedonly if the site to be enteredis empty, exceptfor
external field) is often applied to every particle [25], to the reservoirs i = 0,N +1. After the a step is chosen
approximateaconvectivevelocity. However,particlesdo and particles moved with the assigned probabilities, an-
notexperiencelocal“pondermotive”forces(asisthecase otherparticleispickedatrandomandthepossiblemoves
in osmosisorpressure drivenflows)in the absence ofex- made. ThisalgorithmisidenticaltoaMonteCarlosimu-
ternal electric or gravitational fields. In confined nearly lation[30]exceptthatspatialmovesarecoarse-grainedto
1D systems, such external forces may yield qualitatively ℓ and the effective interaction energies (which determine
different behavior (such as shocks in the particle density the acceptance of moves) between particles are infinite
[28,29])from that expected in hydraulic or osmotic flow. (representing hard spheres) for a particle occupying the
Inthefollowing,weanalyticallymodelparticledynam- neighboringsite,orzeroifthesiteisempty. Thesekinetic
ics through a single pore with the assumption that the steps can be summarized by the following instantaneous
reservoirs are in thermodynamic equilibrium, as was as- fluxes
sumed in the expansion of J (Eqn. (2.1)) about a ther-
modynamicallydefined∆µ0. Weincludesufficientmicro- Ji(t)=pσi(t)(1−σi+1(t))−qσi+1(t)(1−σi(t))
scopic detail within the pore, while carefully obeying all 1<i<N −1
thermodynamic constraints in the connected reservoirs.
J (t)≡J (t)=α(1−σ (t))−γσ (t) i=0
(L)→1 0 1 1
B. One-dimensional Pore Exclusion Model
J (t)≡J (t)=βσ (t)−δ(1−σ (t)) i=N
N→(R) N N N
Our single-file model is similar to but more general (2.2)
than the one vacancy models of Hernandez [12] and
Kohler [13]. The approach is fully microscopic, and Thesedynamicsimplicitly include effectiveexcludedvol-
does not rely on undefined hydrodynamic parameters ume interactions between solvent particles, which are
(on these size scales) such as drag factors, viscosity, etc. strong in confined geometries such as 1D pores. How-
Nonetheless, results from the microscopic theory can be ever, we assume that solvent-solvent attractions within
associatedanddirectlycomparedwiththosefromexper- the pore are weak compared to solvent-pore particle at-
imental measurements. Consider the pore shown in Fig. tractionssuchthatconcertedclusterhopsarerare. These
1whichdepictsasingle-fileornearlysingle-filepore. The processes can be be treated numerically [31] and is out-
maininteractionsthatweincludearetheexcludedparti- side the scope of the present study.
cle(solvent)volumeswithinthepore;thusweadaptone- To extract the steady-state currents and occupations,
dimensionalexclusionmodels[28,16]. Theporeisdivided we take the time average h...i of (2.2) to find
into N sections labelled i, each of length ℓ. The occupa-
tion variable σ (t) defines whether pore section i at time
i
3
hJ (t)i=phσ (t)i−qhσ (t)i+(q−p)hσ (t)σ (t)i J α
i i i+1 i i+1 σ =σ −(i−1) =
i 1
p α+β
hJ (t)i=α(1−hσ (t)i)−γhσ (t)i
0 1 1 p(αβ−γδ)
−
(α+β)[(N −1)(α+γ)(β+δ)+p(α+β+γ+δ)]
hJ (t)i=βhσ (t)i−δ(1−hσ (t)i)
N N N
(i−1)(αβ−γδ)
(2.3) − .
(N −1)(α+γ)(β+δ)+p(α+β+γ+δ)
Exact solutions to the steady state hJ(t)i under cer- (2.7)
tain parameter regimes and in the thermodynamic limit
havebeen foundusing matrixalgebratechniques[28,29]. For an infinitely thin membrane (N = 0), the particles
However, in osmosis and hydraulic pressure driven flows passthemathematicalsurfacebasedupontheirreservoir
across symmetric pores, where the particles do not ex- kinematics alone and never interact with each other or
perience intrinsic or pondermotive forces, and have not the membrane interior while crossing the barrier. Over-
developed collective stream velocities (c.f. Conclusions), taking of particles, when two adjacent, indistinguishable
q = p. That is, a particle is as likely to hop to the particles switch positions, does not contribute to the
left or right if both sites to the left and right are unoc- overall current. Therefore, the model above is also valid
cupied. The only “force” driving osmosis is the direc- as long as each section i is not too wide as to contain
tionallyasymmetricexcludedvolumealongthe lengthof more than one particle at any time. Pore diameters of
thepore,whichisultimatelydeterminedbytheentrance up to 2 ∼3 times the particle diameters will on average
and exit rates at the boundaries i = 1,N. Geometri- preclude multiple occupancy.
cally asymmetric pores would have p 6= q and need to
be treatednumerically. The convectivestreamflow limit
where also p 6=q, is discussed in the Conclusions. Many C. Parameter Relationships
molecular dynamics simulations designed to study flow
in narrow channels impose a microscopic pondermotive The parametersξ ≡(α,β,γ,δ,p)usedinthe aboveki-
force with p 6= q [25] (via a gravitational potential for netic model are conditional probabilities and reveal the
example) in order to accelerate the flow. However, grav- structure and microscopic mechanisms of osmosis and
ity and pressure forces, although macroscopically both pressure driven flows. Here, we explicitly relate ξ to
derived from potentials, do not manifest themselves mi- macroscopicthermodynamic variables. Rather than find
croscopicallyinthesameway. Whenq =pthehigheror- precise numerical values for ξ, we estimate them using
der correlation hσi(t)σi+1(t)i vanishes in (2.3) and upon reasonable physical models and extract their most sen-
time averaging and invoking steady state particle con- sitive dependences, thereby illustrating the qualitative
servation, Ji = Ji+1 = Ji+2 = ... = J0 = JN, where features of osmosis and pressure driven flows through a
we have droppedthe h...i notation. Since all the steady pore.
state currentsacrosseachinteriorsectionisidentical, we The main assumption we make is that the entire sys-
can sum them up along the chain to obtain tem is under local thermodynamic equilibrium (LTE):
p Locally, particle velocities obey a Maxwellian distribu-
J = (σ1−σN)=J0 =JN (2.4) tion. This is justified since measurements of almost all
N −1
microscopic systems give single pore osmotic transport
Equations (2.4) give three equations for the three un- rates of J < 107/s. Typical pore diameters in aqueous
knowns J, σ1, and σN which are found to be (with solutions give mean free paths λMFP . 5˚A. Thus, am-
ξ ≡(α,β,γ,δ,p)) bient thermal velocities v ≃ k T/m ≃ 4×104cm/s
T B
(where m is the particlemass),give meancollisiontimes
p(αβ−γδ) of τ ≃λ /v ≃0.5ps≪pJ−1. Therefore, particles
J(ξ,N)= coll MFP T
(N −1)(α+γ)(β+δ)+p(α+β+γ+δ) within the pore equilibrate by suffering O(105) collisions
while being transported from (L) to (R). This is suffi-
(2.5)
cientforlocalthermodynamicequilibrium(LTE)tohold.
Figures 1(b) are schematics depicting the activa-
and
tion energies E determined by microscopic membrane-
ξ
α−J(ξ,N) δ+J(ξ,N) solvent interactions. Although these are predominantly
σ = , σ = , (2.6)
1 α+γ N β+δ of an enthalpic nature, determined by instantaneous in-
termolecular potentials, a small entropic contribution
which determines σi, mayariseduetoforexamplerotationalaveragingofnon-
spherically symmetric particles. We do not attempt to
model in detail all the configurational dependences of
particletrajectories,butratherassumetheactivationen-
ergiestobeeffectivequantities,averagedoversay,molec-
4
ular rotations, internal vibrations, etc, which will not height E is given by
ξ
qualitativelyaffectthemeanquantitieswewishtostudy.
E areinternalbindingenergiestotheporewallswhich τ ≃[ξ exp(−E /k T)]−1 (2.8)
p,q ξ 0 ξ B
must be overcome when particles hop from one section
to another and can be thought of as effective interac- where ξ0 is a thermal attempt frequency. The transition
tion energies between pore and solvent. As mentioned, probabilities per time ξdt are thus givenby dt/τξ. Upon
the model presented is valid provided solvent-solventat- choosing dt≪τξ,
traction within the pore ≪ E . E and E are the
p β,γ α,δ
effective energies of removing and injecting a particle. If ξ ≃ξ0exp(−Eξ/kBT). (2.9)
we assume that there is no additional “constriction” at
The prefactors of q = p and β,γ, are estimated by
the pore ends, for a pore that repels particles relative to
p ,β ,γ ≃ v /ℓ. Although more sophisticated transi-
itsinteractionenergiesinthebulkreservoirs(L)and(R), 0 0 0 T
tionstatetheoriesexist,thethermalattemptfrequencies
E ≃ 0 and E > 0 (Fig 2(a)). Similarly, as shown
β,γ α,δ are all on the order of 1ps−1. The frequencies α and
in Fig 2(b), an attractive pore will have E > 0 and 0
β,γ
δ however, depend on the thermodynamic state of, and
E ≃ 0. The physical origin of the energies such as 0
α,δ
theparticlenumbersinthereservoirs(L)and(R)respec-
E −E defined in Figure 2(b) are most likely from hy-
α γ
tively, as well as the pore entrance areas. For example,
drogen bonding interactions in aqueous systems. Trans-
for effusing dense gases, α = nLvLAL/4, proportional
ferring a water molecule from the polar, hydrogenbond- 0 T T p
ing bulk environment into the pore region would require to the number density nLT of particles able to enter, the
breaking hydrogen bonds (∼ 5kBT each) and forming thermal velocity vTL, and the available area ALp for en-
van der Waals-type interactions with the pore interior. trance. A completely analogous expression holds for δ0
with (L) replaced by (R). Although the 1/4 factor is
valid for gases, a qualitatively similar expression holds
Eβ =Eγ =0 for dense liquids.
The microscopic pore area A can be defined by
Eα =Eδ L=N l (a). E . First consider E (~r) to bepa general interaction
α,δ α
of pushing a particle perpendicularly through a mem-
brane. When the particle is pushed through at a posi-
tion ~r away from the pore (through a lipid layer for ex-
ample), E (~r) ≃ ∞. When the trajectory runs through
α
(b). the center of an approximately circular pore entrance,
Eα =Eβ =Eγ =Eδ =0 E ≡ E (~r = 0). Therefore, we can define a pore ra-
α α
dius (A ≃ πr2) by demanding activation energies such
p p
that trajectories outside of r are energetically rare. For
p
example,
E =E =0
α δ
Eβ =Eγ (c). Eα(r <rp)−Eα(0)<∼10kBT. (2.10)
Condition (2.10) simply contrains r < r to be within a
p
region where entrance is not impeded by “wall” interac-
FIG. 2. Simplified pore-solvent interaction energy land- tions. Thus, Eα (and Eδ) are the activation energies for
scapefor(a)repulsiveand(b)attractivepores. Theprogres- a particle draggedthrough the approximate “centerline”
sion of the activation energies E are indicated as pores are of the pore. When more than one species occupies the
ξ
altered orconditionsmodifiedsuchthatpore-solventbinding reservoirs,but only one (the solvent) can enter the pore,
increases: (a)→(b)→(c). the entrancerateswillalsobe proportionalto theirmole
fractions. Fornotationalsimplicity,we henceforthdefine
Additional activation energies at the entrance/exit α ≡ α χLexp[−E /k T] and δ ≡ δ χRexp[−E /k T].
0 0 α B 0 0 δ B
siteswouldrepresentforexampleenergeticsofapossible Thus, α ∝ A is now the rate prefactor for pure sol-
0 p
intermediate state where hydrogen bonds are being bro- vententrance (χ =1)at a reference number density for
0
ken, but attractive van der Waals interactions have not a pore with entrance area A . Only the solvent fraction
p
yet formed. For simplicity of notation, we assume these canentertheporeandpartakeinthedynamicsdescribed
are small and neglect them. Inclusion of these bound- by(2.2)becausethesolutesaretoolargetoenterthrough
ary activation barrierswill only shift the absolute values A .
p
of Eα,Eβ,Eγ,Eδ and will not qualitatively change our Although q =p, a net current will occur due to differ-
results. encesin α6=δ, and/orβ 6=γ. Henceforth,for simplicity,
The probabilities ξ can be estimated by first passage weconsideronlystructurallyrightcylindricalporessuch
time or transition state theory. The average time τξ for that β0 = γ0 and α0 = δ0 (or ALp = ARp and Eα = Eδ)
a thermally agitated particle to climb over a barrier of under equilibrium (TL = TR and PL =PR) conditions.
5
Define from the numerator of J(ξ,N) the dimensionless pore is repelling, the energy landscape is approximately
driving variable that shown in Fig. 2(a), where E ≃ E > 0, and
α δ
E ≃ E ≃ 0. As pores with increasing solvent affini-
β γ
αβ−γδ tiesareused,E ≃E decreasewhileα, δ increase,until
X ≡ α δ
αβ the energy barrier shown in Fig. 2(b). is reached, where
α¯ = α /β . Further increasing pore-solvent attraction
χR 0 0
=1− χ0L exp[(Eα+Eβ −Eγ −Eδ)/kBT] (2.11) requires Eα ≃ Eδ ≃ 0, and increasing Eβ ≃ Eγ. Thus,
0 as pores follow the sequence (a) → (b) → (c), α¯ first
=1−e∆E/kBT + ∆χs e∆E/kBT increases as exp(−Eα/kBT) (as Eα > 0 is decreased),
1−χL until α¯ = α0/β0, where upon α¯ increases according to
s
exp(E /k T) (as E > 0 is increased). While the pore
β B β
where∆E ≡Eα−Eγ+Eβ−Eδ ≡ER−EL,∆χs ≡χRs − is repelling, and β =β0 is constant, J(ξ,N) has a maxi-
χL, and χ(R),(L) is the solute mole fractions in (R),(L). mum as a function of pore-solvent affinity at
s s
The solute particles, unlike the solvents, are too large
α∗ 2p/β +(N −1) 1/2
to enter the pore. There will be a steady state particle α¯∗(β )≡ = 0 . (3.2)
0
flux as long as X 6= 0, which can occur either by virtue (cid:18)β0(cid:19) (cid:20) (N −1)(1−X) (cid:21)
of ∆E 6=0 (predominately pressure driven flow), and/or
∆χ ≡χR−χL 6=0 (osmosis). However, if α0 < α∗(β0), the maximum defined by (3.2)
s s s is not reach by decreasing E . Further increasing α/β
Although electrostatic potential differences can con- α
requires decreasing β < β with fixed α = α as shown
tribute to ∆E , they would also make q 6= p, which we 0 0
ξ
bythetransition(b)→(c)inFig. 2. Thevalueofβ which
do not consider. Therefore, we assume that ∆E arise
ξ
yields a maximum in J in this case is
solely to a hydrostatic pressure difference in the baths
which pushes the molecules againsttheir repulsive inter-
1/2
α α (N −1)
actions,andchangetherelativeenergiesofactivatedpore α¯∗(α )≡ 0 = 0
entrance. 0 (cid:18)β∗(cid:19) (cid:20)α0(N −1)(1−X)+p(2−X)(cid:21)
(3.3)
III. RESULTS AND DISCUSSION p=0.003, X=χ=0.0018 [=100mM]
s
12.0
In this sectionwe examine the solutions of one dimen-
(a).
10.0
sional models and discuss their physical meaning. Vari-
ous effects of physical and chemical parameters, such as
pore length and solute interaction dependences, are out- −1) 8.0
s
lined. µ
(
) 6.0
5
1
=
A. Osmosis N N=15 4.0
ξ,
(
J
2.0
Consider symmetric pores (β = γ, E = E ) connect-
α δ
ing two reservoirs under identical hydrostatic pressures
sothat∆E =0. Forsimplicity,we willalsoassumepure 25.0 0.0
−2.0 −1.0 0.0 1.0 2.0
solvent in (L) such that ∆χ ≡ χR. Even though par-
s s
ticle enthalpies are identical for (L) and (R), α0 6= δ0 20.0 (b).
will drive an osmotic flow from (L) → (R) since in this
limit X ≡ (α −δ)/α ≃ χR 6= 0 is due solely to a so-
s 15.0 )
lute concentration difference. The steady state current −1
s
becomes µ
10.0 N=5 ) (
α¯pX N
J(ξ,N)≃ , ξ,
(N −1)(α¯+1)(α¯+1−α¯X)+p¯(2α¯+2−α¯X) (
J
5.0
(3.1)
N=25
where α¯ ≡α/β. Since α/β is a ratio of pore entrance to 0.0
−2.0 −1.0 0.0 1.0 2.0
exit rates across the pore, it measures an effective pore- log (α/β)
solvent binding constant. 10
How does J(ξ,N) behave as pore-solvent combina-
tions with varying α¯ are used? When a symmetric
6
FIG.3. Fluxasafunctionofeffectivepore-solventbinding
0.0 2.0 4.0 6.0 8.0 10.0
α¯. We use representative values dt ≃ 10fs, p = β0 = 0.003, 1.0
and X = χR ≃ 0.0018, corresponding to osmosis from pure
s
solvent to a 100mM solute solution. (a). J(ξ,N = 15) is 0.8 (a).
shown for α0 = 0.001(solid curve) and α0 = 0.0045(dashed
curve) which give difference maxima at α¯∗(α0) and α¯∗(β0) 0.6 σN
respectively. (b) shows J for N = 5,10,15,20,25 when 0.4 , 1
σ
α0 =0.001. Note thediscontinuities in slope at log(1/3) cor- σ(X=0.0018,0.18), σ (X=0.0018)
0.2 1 N
respondingtothepointwherethesourceofincreasingaffinity σ (X=0.18)
switchesfromtheincreasingofαtothedecreasingofβ. (Fig- 0.0 N 0.05
ure2(b).)
(b). 0.04
)
Figure 3 shows the current as the solvent-pore affinity σN 0.03
is increased according to α¯, as α¯ = α /β → α /β (≡ −
0 0 0 1
α¯ ) → α/β . The solute concentration is represented σ 0.02
0 0 (
in molar units [c ] = χR × 55.56, such that a 100mM 0
s s 0 0.01
aqueoussolutesolutioncorrespondstoX=χR ≃0.0018. 1
s
Although J and ξ can all be rescaled by the constant p, 0.00
0.0 2.0 4.0 6.0 8.0 10.0
for concreteness and comparison with experiments, we
α/β [=affinity]
consider a time step dt ≃ 10fs ≪ τ ≈ 0.3ps, which
p
sets β ≈ p ≃ 0.003 and yields flows in the range of
0
108×J (particles per µs) in good agreementwith P ∼ FIG. 4. Occupations σ for p = β0 = α0 = 0.003. (a)
os
10−14cm3/s measured across aquaporin water channels Pore occupations for X =0.0018 are all indistinguishable on
with physiological(∼100mM) osmolyte solutions [2,4]. this scale (solid curve). Only at X = 0.18(10M solution)
Figure 3(a) plots the flux (µs−1) as a function of is σN=30 (dashed curve) sufficiently different from σ1. (b)
logα¯ for two different values, α /β = 0.001/0.003 < 100×(σi−σN) as a function of α/β for p = 0.003,0.03,0.3.
0 0
(α/β )∗ = 8/7(1−X) and α /β = 0.0045/0.003 > Note the different behavior when α = α0 = 0.003 at
0 0 0 α/β = 1.0. The dashed fragments represent the occupation
(α/β )∗. At log(1/3), the curve bifurcates into the two
0 p differences if α0 >0.003.
solutionscorrespondingtoeitherdecreasingαorincreas-
ingβ,withmaximadeterminedby(3.2)and(3.3)respec-
Figure 4(a) shows the average occupation σ (α¯) for
tively. The discontinuity in slope at log(1/3) (dashed i
X = χ = 0.0018. The dependence on spatial position i
curve) is apparent, although the one at log(3/2) near s
is weakandis notapparentonthe scaleofFig. 4(a), ex-
the maximum is not. Figure 3(b) compares J(ξ,N) for
cept if X = 0.18 (corresponding to 10M solutions) when
various pore lengths L = Nℓ at α = 0.001 (such that
0 σ can be distinguished from σ . Therefore, from Eqn.
(3.3) pertains). Whether α /β is greater or less than N 1
0 0 (2.6) α/(α+β) is an approximate value of averaged oc-
α¯∗(β ) depends onthe reservoirnumber density and A ;
0 p cupations at small X. The small difference in occupa-
for gases, α ≪ β and a curve with maximum defined
0 0 tion along the chain is consistent with the assumption
by (3.3) obtains. However, dense liquids entering large
of local thermodynamic equilibrium. Figure 4 shows the
pores may have relatively large entrance rates such that
difference 103×(σ (α¯)−σ (α¯)) with α /β =1.0(solid
α /β & α¯∗(β ) and J(ξ,N) may have a maximum de- 1 N 0 0
0 0 0 curves). Thedashedcurvesrepresenttheoccupationdif-
fined by (3.2). From (3.2), this latter condition is less
ference 103×(σ −σ ) if α /β ≫1.
likely for larger driving X, and thus α¯∗(β ). 1 N 0 0
0 The behaviorof J(ξ,N) asa function of affinityα¯ can
be further understood in terms of the associated pore
occupations. At low α¯, the pore is repelling and con-
tains few if any particles that can be transported. At
high affinities, and occupations, most internal hopping
steps are not fulfilled due to occupied neighboring sites
choking off the flow. Only at intermediate α¯ and near
(3.2) or (3.3), and intermediate occupations will maxi-
mum flow occur. At larger X, σ decreases (see Fig.
N
4(a)), decreasing overall pore occupation and increasing
the critical affinities α¯∗(β ) and α¯∗(α ) where the max-
0 0
imum occurs. Therefore, the condition represented by
(3.3)becomesmorerelevantforlargerosmoticpressures.
7
dt=10fs, p=β =0.003, α =0.002
From the denominator of (3.1), J is linear in X for
0 0
25 500 α¯X/(α¯ +1) ≪ 1. From (2.6), we see that α¯/(α¯ +1) is
(a). (b). the average occupation for small J. Therefore, nonlin-
450
earities are important only when occupation σ is high,
i
20 400 when exclusion interactions, and thus nonidealities, are
most pronounced (the α¯ = 1.0,100 curves in Fig. 5(b)).
350
Although Eqn. (3.1), or (3.4) fully describes particle ex-
15 1) 300 clusion nonidealities within the pore, for the small ∆χs
−s usually encountered, a linear relationship,
µ
( 250
) α¯ p∆χ
0 s
10 2 200 J(ξ,N)≃ (3.5)
ξ, α¯+1(N −1)(α¯+1)+2p¯
(
J 150
is accurate as demonstrated by Fig 5(a). In the limit
5 100 p¯≪(α¯+1)N, Pos approaches
50 P ∼= a αp¯ . (3.6)
os L (α¯+1)2n
T
0 0 (cid:16) (cid:17)
0.0 0.2 0.4 0.6 0.8 1.0 0 5 10 15 20
In this limit(3.6), the “rate limiting steps” are the
X=χ (Molar)
internal steps (p) of pore particles, and the overall
s
flux scales as L−1. In the limit represented by (3.6),
FIG.5. ThedrivingforcedependenceofJ(ξ,20)forosmo-
an Arrhenius temperature dependence of (E − E −
sis from pure solvent to solution with p = β0 = 0.003 and α β
E )/k T [(E −E −E )/k T] is expected for α¯ ≫
α0 = 0.002. Various values of α¯ = α/β are shown. (a) For q B β α q B
solute concentrations .5M, the current is linear in X=χR. 1[α¯ ≪1].
s Accordingtotheformofα ,P ∝αp¯∝r2 forα¯ ≪1,
(b)Athighersoluteconcentrations, nonlinearities inJ(ξ,20) 0 os p
arise, particularly for high valuesof α¯ an pore occupation. and Pos ∝ p/α¯ ∝ rp−2 for α¯ ≫ 1. Note the curious
P ∝ r−2 dependence indicating the α¯ > α¯∗ regime
os p
Figures 5 show the particle flux per microsecond, shown by the large log α¯ regions of Fig. 3 where de-
10
108 × J(ξ,N = 20), as a function of X ≃ χR at creasing r (and α) increases J¯ by decreasing the high
s p
β = p = 0.003, α = 0.002 for various pore-solvent pore occupation.
0 0
affinities α¯. For the low values α¯ = 0.01(dashed curve) Nowconsiderthelimitp¯≫(α¯+1)N. Theratelimiting
andα¯ =0.1(dottedcurve)theporeisrepelling,whilefor steps now involve pore entrance or exit and
α¯ = 1.0(solid curve) and α¯ = 100(thick curve) α = α
0 α
while β decreases from β0. The solid curve with α¯ =1.0 Pos ∼= , (3.7)
2(α¯+1)n
is near the maximum defined by (3.3) while α¯ = 100 is T
within the choked flow limit represented by the decreas-
independent of L. This limit is equivalent to the single
ingbranchonthehighα¯ sideofFigures3(a,b). Theflow
site(N =1)modelsinceonlytheboundaryentranceand
is linear for low X as shown in 5(a). Figure 5(b) clearly
exit rates are relevant. The radius and temperature de-
shows that nonlinearities are important for large α¯ and pendences in this case are P ∝ α ∝ r2exp(−E /k T)
that the flow for α¯ = 1.0(solid curve) is lower that that os p α B
and P ∝ β ∝ exp(−E /k T) for α¯ ≪ 1 and α¯ ≫ 1
os β B
for α¯ = 0.01(dashed curve) only for small driving X. As
respectively. Thus, the pore length L ≃ Nℓ can deter-
the nonlinearity sets in for α¯ = 1.0, the flux increases
mine the temperature dependence of J since it deter-
fasterthanthatforα¯ =0.01andsurpassesitnearX∼3.
mines which of the two limits (3.6) or (3.7) is relevant.
The flow nonlinearities shown in Fig 5(b) arise at large
Exceptforexperimentswherelipids undergophasetran-
soluteconcentrationdifferencesandareonlythoseresult-
sitionswhichaffecttheporestructure[33],asingleexpo-
ingfromthenonlineardynamicsdeterminedbyequations
nentialtemperature dependence is almostalwaysexperi-
(2.2) and do not include solute nonidealities, unstirred
mentallyobservedinbiologicalsystemsinvolvingosmosis
layers, etc. That is, the coefficients
[32], indicating one of the above limits for P holds in
os
α¯kp[(α¯+1)(N −1)+p¯]k−1 the regimes explored. One interpretation is that a num-
dk = (α¯+1)k[(α¯+1)(N −1)+2p¯]k (3.4) ber of hydrogen bonds (∼ 5kBT each) must be broken
before water can enter the pore. This is consistent with
contain only the kinetic parameters describing the the identification of E as the energy required to break
α
solvent-pore interactions. Furthermore, our treatment waterH-bondsintheα¯ ≪1limitof(3.6). Fluxesofnon-
for α and β assume thermodynamic equilibrium, or no polar solvents through artificial pores may be expected
unstirred/polarization layers. The modifications neces- to display richer temperature dependences since a wider
saryto include unstirredlayersaretreated inthe follow- variety of temperatures and molecular interactions can
ing paper. be accessed and experimentally probed.
8
B. Pressure Driven Flow where in this case,
X=1−exp (E (PL)+E −E (PL)−E )/k T ,
For incompressible reservoirs without solutes, ∆χ = α β γ δ B
s
χLs =χRs =0; however,flow can be driven by differences (cid:2) (3(cid:3).11)
in hydrostatic pressure. Hydrostatic pressure variations,
rather than controlling the rates via the number density andd¯measuresthesolvent-poreaffinitysinceα¯isvarying
in the pre-exponentialfactors α0, δ0, modify the relative with pressure. For low pressure differences such that
activation energies E . For example, consider increas-
ξ
ing PL keeping PR constant. For nearly incompressible −1
∂∆E
liquids, the particle density changes negligibly, but their PL−PR ≪Eδ (3.12)
∂∆P
particle enthalpies increase as they are pushed against "(cid:18) (cid:19)T,∆P=0#
their interaction potentials. Changing the relative reser-
flowacrossrepellingporescanbefurthersimplifiedsince
voir and pore enthalpies changes the entrance and exit
γ =β as evident in the lower three pressures PL in Fig.
kineticcoefficientsaccordingtoFig. 6(a),(b)forinitially
6(b). When the pressure pushes the enthalpy of the left
repelling or attractive pores.
bathabovethatofthepore,thefullpressuredependence
(3.10) is required. Equations (3.8) and (3.10), along
PL E α =0 with the equation of state of the solvent in the reser-
Eδ =0 voirs completely determine the nonlinear flux-pressure
E γ (P L )>0
Eβ >0 (a). relationship. However, for low pressures, a Taylor ex-
pansion about PL =PR yields
∆P ∂∆E ∆P
X≈− =−v˜ , (3.13)
k T ∂∆P k T
PL Eα =0, Eγ (P L )>0 B (cid:18) (cid:19)T,∆P=0 B
E β =0, E δ >0 where the lastequality is valid for incompressibleliquids
(b). andv˜isthemolecularvolumeofasolventparticleinthe
E γ =0, Eα (P L )>0 reservoirs(forH2O,v˜≃3×10−23cm3). Linearizing(3.8)
and(3.10)aboutsmallpressuredifferences,andrecalling
FIG.6. Hydrostaticpressureincreasesinleftreservoir. (a) that α≃δ, we find
In attractive pore, α=δ, Eβ is constant, and Eγ =Eγ(PL).
(b) For repelling pore, and small pressure differences (lower P
solidanddottedlevels)γ =β, Eδisfixed,andEα=Eα(PL). Lp ∼= k oTs (3.14)
Pressure differences larger than Eδ result in Eα ≈ 0 and B
Eγ =Eγ(PL)>0. as required by linear nonequilibrium thermodynamics.
However, differences between hydraulically driven and
Since we assume only PL is varied, E , E are keptcon- osmotically driven transport arise at higher order in X,
δ β
stant, while E = E (PL) and E = E (PL). First as explicitly calculated in (2.5), (3.8), and (3.10).
α α γ γ
consider an attractive pore where only E (PL) and γ Continuum Poiseuille flow (zero Reynolds number so-
γ
changes as PL is increased, as schematically depicted in lution of the Navier-Stokes equation for flow through an
Fig. 6(a). The flux is infinite right circular cylinder [19]) expressions (where
J ∝ r4) have been proposed to describe transport
p
J γ(PL),N = through microscopic channels. Although these contin-
pα¯X (3.8) uum relations have no validity in microscopic settings,
(cid:0) (cid:1)
(N −1)(α¯+1−X)(α¯+1)+p¯(2α¯+2−X) they are nevertheless often used to obtain estimates of
permeabilities, pore radii [20], and membrane pore den-
where sities. These estimates typically do not agree well with
independentmeasurementsofporeradiiandporeperme-
X(PL)≡1−exp[∆E/k T] (3.9)
B abilities [20]. Although the r dependences in our model
p
are valid only for a limited range of pore radii, they are
and ∆E ≡E −E (PL).
β γ consistentwith the inapplicability, atmicroscopicscales,
When hydraulic pressure pushes particles through a
of the continuum expression P ∝ r4/(ηL), where η is
repelling pore, the current can be written in the form os p
the bulk solvent viscosity.
J α(PL),γ(PL),N = Fluid flow through small orifices has been described
pδ¯X as being both “viscous” and “diffusive.” Even though
(cid:0) (cid:1) the limit (3.6) appears to be “viscous” due to the L−1
(N −1)(δ¯+1)(δ¯+1−X)+p¯ γ¯+δ¯+(δ¯+1)(1−X)/γ¯
dependence, and and (3.7) “diffusive,” we see that both
(cid:0) (3.10) (cid:1) pressure and osmosis driven flows in the present model
9
is diffusive in nature. Consider the time averaged rate by pσ (1−σ )(1−σ )+p′σ (1−σ )σ where p
i i+1 i−1 i i+1 i−1
equation for occupation at site i, represents the rate of an isolated particle, while p′ ≪ p
represents the rate in the presence of another particle at
σ˙i =−pσi(1−σi+1)−pσi(1−σi−1) i−1 binding to and preventing the hop of the particle
+p(1−σi)(σi+1+σi−1) at i to i+1. The resulting particle density profile along
(3.15)
the chain will no longer be linear. These effects offer
=p(σi+1−2σi+σi−1) the possibility of interesting nonlinearities amenable to
numerical study.
whichonscales≫ℓisadiffusionequationwithpℓ2 ∼=Dc For pores that are slightly larger, with radii larger
a cooperative diffusion coefficient. When particles in- than a few times the particle diameters, similar stochas-
teract along the chain the overall flux is determined by tictreatmentscanbeappliedprovidedparticlemomenta
the cooperative diffusion of particles. In the simple one- relaxmuchfasterthanparticlepositionaldegreesoffree-
dimensionalmodel presented,cooperativediffusionis re- dom. Onlythenisthecollectivestreamvelocityzero,and
lated directly to the precisely defined dynamical rules p=q can be invoked. Near the single file limit, the pore
(2.2), but in general is defined by linear density-density walls are in constant interaction with the transported
correlation functions, or Green-Kubo relations. molecules and randomizing their momentum. Neverthe-
The neglect of true convective,or momentum transfer less, for larger pores, we can estimate the relative rates
effects can be neglected in the systems we consider. The of positional and momentum relaxation. Pressure (ei-
signature of convective flow, such as Poiseuille flow, is a therhydrostaticorosmotic)drivenflowthroughamacro-
particlevelocitydistributionfunctionwithisMaxwellian scopic circular pipe is described by Poiseuille’s law [19]
centered about a small finite stream velocity. However,
in 1D, or nearly 1D pores, the particles thermally equi- v(r)= ∆P (r2−r2) (4.1)
librate with the stationary pore walls very quickly and 4ηL p
are unable to collectively transfer a stream momentum.
where v(r) is the fluid element velocity along the pipe,
This is particularly valid for repelling pores where en-
and η is the bulk fluid viscosity. Equation (4.1) has of-
tropic pressure results in low pore occupations σ ≪ 1.
i ten been extended to describe flow through microscopic
In this large Knudsen number regime [6], the particles
pipes,especiallybiologicalmembranepores,havingradii
will collide with walls more often that with neighboring
r of Angstroms. This is the origin of the notion that
particles. p
osmosisthroughmicroscopicporesis“convective.” How-
ever,forthe Poiseuillestreamconvectiontobecompara-
ble to “driven diffusive transport,” the ratio
IV. EXTENSIONS AND CONCLUSIONS
vL ∆Pr2
An extension of the 1D chain represented by (2.2) to Pe≡ ≃ p ∼1, (4.2)
D 4ηD
c c
include a distribution of site-dependent internalhopping
ratespi =qi canbe straightforwardlycalculated,as out- where Dc is the effective collective diffusion coefficient
lined in the Appendix. Biological realizations of inter- which describes dissipative particle density relaxation,
nal pore defects with site dependent hopping rates are and is also approximately the momentum diffusion, or
gramicidin A channels comprising of two barrel struc- soundattenuationcoefficient. For a hydrostaticpressure
tures joined at the defect (membrane bilayer midplane) difference of ∆P = 2.5 atm (corresponding to a 100mM
[1]. Similarly, water permeation through a bilayer mem- solute concentration difference), and bulk viscosity and
brane can be viewed as pores along the hydrophobic attenuation coefficient η = 0.01g cm/s, and Dc ≃ η/ρ,
lipid tails with the bilayer midplane and possible un- convective, or “viscous” [14] flow will be irrelevant for
saturated bonds acting like defects. Channel proteins rp ≪ 1000˚A. For typical zeolites and biological pores
with widely varying internal interactions would also be undertheaboveconditions,r ∼3−8˚A,Pe∼10−5 ≪1.
p
modifiedaccordingtotheresultsfoundintheAppendix. Thus, “convective” flow described by (4.1) is negligible
Strongparticle-poreinteractionswhichareincommensu- comparedtodiffusivetransportinsystemsofpresentcon-
rate with ℓ may also lead to nonlinear transport. Here, cern. Diffusive transport in this context should not be
possible cluster diffusion [31] and locked-sliding phases confused with tracer diffusion or collective diffusion in
such as that found in Frenkel-Kontorowa [34] and fric- binary mixtures, rather, it is the consequence of micro-
tion models [35] can also affect transport. These effects scopic momentum relaxation due to frequent collisions
offerthepossibilityofinterestingnonlinearitiesamenable with the pore interior. These collisions increase relative
tonumericalstudy. Solventsthathavestrongmutualat- to particle-particle interactions for smaller and smaller
tractiveinteractionscanalsobeconsidered. Forexample, pores, where the inner surface to pore volume ratio in-
aparticlemayhopwithdifferentratesintoanemptysite creases.
at i+1 depending on whether or not there is an attrac- Additional support for the basic assumption of diffu-
tive particle behind it, at site i−1. The probability of siveparticlemotionsisfoundfromnonequilibriummolec-
suchahoptotherightfromi→i+1canbe represented ulardynamicsstudies(NEMD)offlowthroughaslit[14]
10