Polymer translocation through a nanopore: A showcase of anomalous diffusion
J. L. A. Dubbeldam,1,2A. Milchev,1,3V. G. Rostiashvili,1and T. A. Vilgis1 1
Max Planck Institute for Polymer Research, 10 Ackermannweg, 55128 Mainz, Germany
2
Delft University of Technology, 2628CD Delft, The Netherlands
3
Institute for Physical Chemistry, Bulgarian Academy of Science, 1113 Sofia, Bulgaria
共Received 23 January 2007; revised manuscript received 21 June 2007; published 13 July 2007兲 The translocation dynamics of a polymer chain through a nanopore in the absence of an external driving force is analyzed by means of scaling arguments, fractional calculus, and computer simulations. The problem at hand is mapped on a one-dimensional anomalous diffusion process in terms of the reaction coordinate s共i.e., the translocated number of segments at time t兲 and shown to be governed by a universal exponent␣=2/共2 + 2 −␥1兲, where is the Flory exponent and␥1is the surface exponent. Remarkably, it turns out that the value of␣ is nearly the same in two and three dimensions. The process is described by a fractional diffusion equation which is solved exactly in the interval 0⬍s⬍N with appropriate boundary and initial conditions. The solution gives the probability distribution of translocation times as well as the variation with time of the statistical moments具s共t兲典 and 具s2共t兲典−具s共t兲典2, which provide a full description of the diffusion process. The comparison of the analytic results with data derived from extensive Monte Carlo simulations reveals very good agreement and proves that the diffusion dynamics of unbiased translocation through a nanopore is anomalous in its nature.
DOI:10.1103/PhysRevE.76.010801 PACS number共s兲: 82.35.Lr, 87.15.Aa
The dynamics of polymer translocation through a pore has recently received a lot of attention and appears highly rel-evant in both chemical and biological processes关1兴. The the-oretical cosideration is usually based on the assumption 关2–4兴 that the problem can be mapped onto a one-dimensional diffusion process. The so-called translocation coordinate 共i.e., reaction coordinate s兲 is considered as the only relevant dynamic variable. The whole polymer chain of length N is assumed to be in equilibrium with a correspond-ing free energy F共s兲 of an entropic nature. The one-dimensional 共1D兲 dynamics of the translocation coordinate then follows the conventional Brownian motion, and the one-dimensional Smoluchowski equation 关5兴 can be used with the free energyF共s兲 playing the role of an external potential. In the absence of external driving force共unbiased transloca-tion兲, the corresponding average first-passage time follows the law共N兲⬀a2N2/ D, where a is the length of a polymer Kuhn segment and D stands for the proper diffusion coeffi-cient. The question of the choice of the proper diffusion co-efficient D, and the nature of the diffusion process, is con-troversial. Some authors 关2,3兴 adopt D⬀N−1, as for Rouse diffusion, which yields ⬀N3 as for polymer reptation 关8兴, albeit the short pore constraint is less severe than that for a tube of length N. In Ref.关4兴 it is assumed that D is not the diffusion coefficient of the whole chain but rather that of the monomer just passing through the pore. The unbiased trans-location time is then predicted to vary as⬀N2. The latter assumption has been questioned关6,7兴 too. Indeed, on the one hand, the mean translocation time scales关4兴 as⬃N2, but on the other hand the characteristic Rouse time共i.e., the time it takes for a free polymer to diffuse a distance of the order of its gyration radius兲 scales as Rouse⬀N2+1, where the Flory exponent = 0.588 at d = 3, and = 0.75 at d = 2 关8兴. Thus
RouseⰇ, against common sense, given that the unimpeded motion should be in any case faster than that of a constrained chain. Moreover, the equilibration of the chain is question-able when the expression forF共s兲 is to be used. The charac-teristic equilibration time scales again as eq⬀N2+1 and is
thus always larger than the translocation time, i.e., eqⰇ. Again the internal consistency of the whole approach is in doubt. It was found by Monte Carlo 共MC兲 simulation 关6,7兴 that⬀N2.5for translocations in d = 2. This indicates that the translocation time scales in the same manner as the Rouse time, albeit with a larger prefactor that depends on the size of the nanopore. Kantor and Kardar argued that this finding bears witness to the failure of the Brownian nature of the translocation dynamics and suggested instead that anoma-lous diffusion dynamics 关9兴 should be more adequate. The
⬀N2.5 scaling law has been corroborated by a further MC study 关10兴, as well as by MC simulations on a 3D lattice 关11兴, and it was shown that⬀N2.46±0.03. The time variation of the second statistical moment,具s2典−具s典2⬀t␣, clearly indi-cates an anomalous nature关11兴, since the measured exponent
␣= 0.81± 0.01, while⬀N2/␣. Still missing is a proper theo-retical analysis which could explain the physical origin of the anomalous dynamics, and make it possible to solve the ap-propriate fractional diffusion equation共or, in case of a biased translocation, the fractional Fokker-Planck-Smoluchowski equation兲 关9,12兴 governing this dynamics.
In this Rapid Communication, we suggest a unique physi-cal picture that justifies the mapping of the 3D problem on a 1D reaction cooordinate s, and we show that the latter obeys anomalous diffusion dynamics, described by a fractional dif-fusion equation. We solve this equation exactly, subject to the proper boundary conditions, and find a perfect agreement with our scaling prediction. Eventually, we demonstrate that the results of our 3D off-lattice MC simulations are in accord with our analytical findings.
Mapping onto 1D dynamics. As indicated above, the as-sumption that the whole polymer chain is in equilibrium and the diffusion is governed by conventional Brownian dynam-ics leads to contradictions. Instead, we assume now that only a part of the whole chain may equilibrate between two suc-cessive threadings. This part of the chain which adjoins the membrane on the cis or trans side will be denoted as fold, and we assume that it is much shorter than the whole chain length N but is still long enough so that one can use the PHYSICAL REVIEW E 76, 010801共R兲 共2007兲
RAPID COMMUNICATIONS
principles of statistical physics. We also assume that the ex-cluded volume interaction of a fold with the rest of the chain is relatively weak, so that it could be treated as a subsystem with a well-defined free energy. This latter assumption is based on the observation that the chain on either of the two sides may be seen as a polymeric “mushroom” whereby the monomer density close to the membrane共or wall兲 is much smaller than the density inside a single coil关see Fig. 4 and Eq. 共II.4兲 in Ref. 关13兴兴. Thus one can claim that there is a depletion area near the membrane关14兴.
Figure1illustrates how a fold squeezes from the cis to the trans side through a short nanopore共of length ⬇a兲, which is slightly wider than the chain itself. It is self-evident that, in the absence of the external force, with the equal probability folds from the trans side could go to the cis side. If the trans part of the fold in Fig.1has length n then the corresponding free energy function Ft共n兲/T=−n ln−共␥
1− 1兲ln n, where is the connective constant and␥1is the surface entropic ex-ponent 关15兴. For the cis part of the fold, one has Fc共n兲/T = −共s−n兲ln−共␥1− 1兲ln共s−n兲 so that the total free energy is
F共n兲/T=−s ln−共␥1− 1兲ln关n共s−n兲兴. One can, therefore, as-cribe to the fold cis-trans transition a pretty broad barrier given by F共n兲. The corresponding activation energy of the fold can be calculated as ⌬E共s兲=F共s/2兲−F共1兲=共1 −␥1兲T ln s.
How can we estimate the characteristic time for the fold cis-trans translocation? In the absence of a separating mem-brane this would be the pure Rouse time tR⬀s2+1 关8兴. The membrane with a nanopore imposes an additional entropic activation barrier ⌬E共s兲, which slows down the transition rate. The characteristic time, therefore, scales as t共s兲 = tR共s兲exp关⌬E共s兲兴⬀s2+2−␥1. This makes it possible to esti-mate the mean-squared displacement of the s coordinate:
具s2典 ⬀ t2/共2+2−␥1兲. 共1兲 Hence, the mapping onto the s coordinate leads to an anoma-lous diffusion law具s2典⬀t␣, where␣= 2 /共2+ 2 −␥
1兲. Taking into account the most accurate values of the exponents for d = 3,= 0.588, and ␥1= 0.680 关16兴, we obtain␣= 0.801. In turn, the average translocation time⬀N2/␣⬀N2.496. Remark-ably, in 2D, where2D= 0.75 and ␥1⬇0.945 关17兴, one finds
␣⬇0.783 共i.e., ␣ is almost unchanged兲. This explains why the measured exponents in both 2D关6兴 and 3D 关11兴 are so close. The derivation of␣is our central scaling prediction. It also agrees well共see below兲 with our own MC data on the translocation exponent.
Fractional diffusion equation. We now turn to the frac-tional diffusion equation 共FDE兲 that furnishes a natural
framework for the study of anomalous diffusion关9,12兴. Here we make use of this method in a systematic way. Our FDE reads tW共s,t兲 =0Dt 1−␣K ␣ 2 s2W共s,t兲, 共2兲 where W共s,t兲 is the probability distribution function 共PDF兲 for having a segment s at time t in the pore, and the fractional Riemann-Liouville operator 0Dt1−␣W共s,t兲 =关1/⌫共␣兲兴共/t兲兰0tdt
⬘
W共s,t⬘
兲/共t−t⬘
兲1−␣. In Eq. 共2兲 ⌫共␣兲 is the Gamma function, and K␣ is the so-called generalized diffusion constant. This constant could be defined as K␣ =⌫共1+␣兲l2/共2w
␣兲 in terms of the fold length l and the wait-ing time scalew共see Chapter 3.4 in 关9兴兲. It should be men-tioned that the constant K␣ is the only adjustable parameter of our theory, and will be fixed below through the compari-son with our MC data.
Recently the method of the generalized Langevin equa-tion共GLE兲 has been used to describe anomalous conforma-tional dynamics within single-molecule proteins关18兴. In con-trast to the FDE approach, which deals with the total distribution function at particular boundary conditions 共see below兲, the GLE method treats only the first two moments 共or time-correlation functions, memory kernel, etc.兲. To the best of our knowledge, at the present time it is not clear how one can derive in a closed form a non-Markovian Fokker-Planck equation for the distribution function 关19兴 starting from the GLE. On the other hand, the translocation time distribution function共see below兲 is an entity of great impor-tance because it could be directly measured in experiment 关1兴. Therefore, we prefer to use the FDE approach for the translocation problem.
Consider the boundary value problem for the FDE in the interval 0艋s艋N. This problem has been discussed before in the context of the even more general fractional Fokker-Planck equation关20兴. The boundary conditions correspond to the reflecting-adsorbing case, i.e., 兩共/s兲W共s,t兲兩s=0= 0 and W共s=N,t兲=0. The initial distribution is concentrated in s0, i.e., W共s,t=0兲=␦共s−s0兲. The full solution can be represented as a sum over all eigenfunctions n共s兲, i.e., W共s,t兲
=兺n=0⬁ Tn共t兲n共s兲, where n共s兲 obey the equations
K␣共d2/ ds2兲
n共s兲+n,␣n共s兲=0, and the eigenvalues n,␣can
be readily found from the foregoing boundary conditions; as a result n,␣=共2n+1兲22K
␣/共4N2兲. The temporal part Tn共t兲
obeys the equation 共d/dt兲Tn共t兲=−n,␣ 0Dt1−␣Tn共t兲. The
solu-tion of this equasolu-tion is given by Tn共t兲=Tn共t=0兲E␣共−n,␣t␣兲
关9兴, where the Mittag-Leffler function E␣共x兲 is defined by the series expansion E␣共x兲=兺n=0⬁ xk/⌫共1+␣k兲. At ␣= 1 it turns
back into a standard exponential function共normal diffusion兲. Thus we arrive at the complete solution of Eq.共2兲:
W共s,t兲 =2 N
兺
n=0 ⬁ cos冉
共2n + 1兲s0 2N冊
cos冉
共2n + 1兲s 2N冊
⫻ E␣冉
−共2n + 1兲 22 4N2 K␣t ␣冊
. 共3兲First-passage time distribution. The distribution of trans-location times 共which could, in principle, be measured in
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 (a) (b) s s−n n CIS TRANS TRANS CIS
FIG. 1.共Color online兲 Schematic representation of a chain fold of length s moving through a nanopore. The transition rate is slowed down by an entropic barrier:共a兲 initially the fold is on the
cis side of the wall;共b兲 the fold entropy decreases during threading
because of the fold fragmentation.
DUBBELDAM et al. PHYSICAL REVIEW E 76, 010801共R兲 共2007兲
RAPID COMMUNICATIONS
experiment兲 is nothing but the first-passage time distribution 共FPTD兲 Q共s0, t兲, where s0 stands for the initial value of the
s coordinate. Knowing the probability distribution function W共s,t兲, we can calculate the FPTD Q共s0, t兲. The relation between the two functions is given as Q共s0, t兲 = −共d/dt兲兰0
t
W共s,t兲ds 关5兴. This yields the FPTD as follows:
Q共s0,t兲 =K␣t ␣−1 N2 n=0
兺
⬁ 共− 1兲n共2n + 1兲cos冉
共2n + 1兲s0 2N冊
⫻ E␣,␣冉
−共2n + 1兲 22 4N2 K␣t ␣冊
, 共4兲where the generalized Mittag-Leffler function E␣,␣共x兲 =兺k=0⬁ xk/⌫共␣+ k␣兲. The long time limits of Mittag-Leffler functions in Eqs. 共3兲 and 共4兲 follow an inverse power law behavior, E␣共−n,␣t␣兲⬀1/⌫共1−␣兲n,␣t␣ and E␣,␣共−n,␣t␣兲
⬀␣/⌫共1−␣兲n,2␣t2␣. By making use of this in Eq. 共4兲, the long time tail of the FPTD then reads Q共t兲⬀␣N2/ 2⌫共1 −␣兲K␣t1+␣. This behavior is checked below in our MC in-vestigation. It can be seen that the mean first-passage time, defined simply as =兰0⬁tQ共t兲dt, does not exist 关21,22兴. On the other hand, in a laboratory experiment there always ex-ists some upper time limit t*. Taking this into account, one can show that an “experimental” first-passage time scales as
⬃N2/␣关21兴, which we observe in our MC simulation.
Statistical moments具s典 and 具s2典 vs time. The subdiffusive behavior of the second moment具s2典−具s典2⬀t␣ is a hallmark of anomalous diffusion. Starting from Eq.共3兲 we can imme-diately calculate them. The calculation of the first moment 具s典=兰0 N sW共s,t兲ds/兰0NW共s,t兲ds yields 具s典共t兲 N = 1 − 2
兺
n=0 ⬁ 1 共2n + 1兲2E␣冉
− 共2n + 1兲22 4N2 K␣t ␣冊
兺
n=0 ⬁ 共− 1兲n 共2n + 1兲E␣冉
− 共2n + 1兲22 4N2 K␣t ␣冊
. 共5兲 Since E␣共t=0兲=1, the initial value 具s典共t=0兲=0 共we put s0 = 0兲, as it should be. In the opposite limit, t→⬁, we can use the asymptotic behavior E␣关−n,␣t␣兴⯝1/⌫共1−␣兲n,␣t␣as well as the sum values 兺n=0⬁ 1 /共2n+1兲4=4/ 96 and 兺n=0⬁ 共−1兲
n/共2n+1兲3=3/ 32 in the numerator and denomina-tor, respectively. After that具s典共t→⬁兲=N/3, i.e., the function goes to a plateau.
The result for the second moment 具s2典 =兰0Ns2W共s,t兲ds/兰0NW共s,t兲ds can be cast in the following form: 具s2典共t兲 N2 = 1 − 8
兺
n=0 ⬁ 共− 1兲n 共2n + 1兲3E␣冉
− 共2n + 1兲22 4N2 K␣t ␣冊
2兺
n=0 ⬁ 共− 1兲n 共2n + 1兲E␣冉
− 共2n + 1兲22 4N2 K␣t ␣冊
. 共6兲Again, it can be readily shown that 具s2典共0兲−具s典2共0兲=0. In the long time limit, in the same way as before and taking into account that 兺n=0⬁ 共−1兲n/共2n+1兲5= 55/ 1536, we find 具s2典共t→⬁兲−具s典2共t→⬁兲=N2/ 9 关23兴.
Monte Carlo data vs theory. We have carried out exten-sive MC simulations in order to check the main predictions
of the foregoing analytical theory. We use a dynamic bead-spring model which has been described before关24兴; therefore we mention only the salient features here. Each chain con-tains N effective monomers 共beads兲, connected by anhar-monic finitely extensible nonlinear elastic springs, and the nonbonded segments interact by a Morse potential. An el-ementary MC move is performed by picking an effective monomer at random and trying to displace it from its posi-tion to a new one chosen at random. These trial moves are accepted as new configurations if they pass the standard Me-tropolis acceptance test. It is well established that such a MC algorithm, based on local moves, realizes Rouse model dy-namics for the polymer chain. In the course of the simulation we perform successive run for chain lengths N = 16, 32, 64, 128, 256, whereby a run starts with a configura-tion with only a few segmens on the trans side. Each run is stopped, once the entire chain moves to the trans side. Com-plete retracting of the chain back to the cis side is prohibited. During each run we record the translocation time and the translocation coordinate s共t兲. Then we average all data over typically 104 runs. In principle, the pore may apply a drag force on the threading chain due to a chemical potential gra-dient; however, in the present work we consider only unbi-ased diffusion. In Fig. 2共a兲 we show a master plot of the translocation time distribution Q共t兲 derived from Eq. 共4兲 for different chain lengths N = 16, 32, 64, 128, 256. For the calcu-lation of data we have used MATHEMATICA with a special
package for computation of Mittag-Leffler functions 关25兴. Evidently, all curves collapse on a single one when time is
10−2 100 t/N2ν+2−γ1[MCS] 10−4 Q(t) N=16 N=32 N=64 N=128 N=256 t−1.8 0 25000 50000 t [MCS] 0 2e−05 4e−05 6e−05 Q(t) (a) 101 102 103 t/N2ν+2−γ1 10−3 10−2 Q (t) N=64 N=128 τ−1.8 101 102 N 104 106 108 <τ (N)> N2.52 +/− 0.04 (b)
FIG. 2. 共Color online兲 Translocation time distribution function
Q共t兲. 共a兲 Scaling plot of the theoretical predictions calculated from
Eq.共4兲 for different chain lengths N. Dashed line denotes the long
time asymptotic tail with slope −1.8. The inset shows Q共t兲 for N = 256 in normal coordinates.共b兲 The FPDT Q共t兲 from the MC simu-lation for N = 64, 128. The inset shows the expected具典 vs chain length N dependence, and the straight line is a best fit with slope ⬇2.52±0.04.
POLYMER TRANSLOCATION THROUGH A NANOPORE: A… PHYSICAL REVIEW E 76, 010801共R兲 共2007兲
RAPID COMMUNICATIONS
scaled as t⬀N2/␣ with the predicted ␣= 0.8. The long time tail for this value of ␣ should exhibit a slope of −1.8. The inset in Fig.2共a兲reveals the long tail of Q共t兲 for large times. A comparison with Fig. 2共b兲 demonstrates good agreement
with the simulation data despite some scatter in the FPDF even after averaging over 10 000 runs. As shown in the inset, the mean translocation time scales as 具典⬀N2.5 in good agreement with the predicted␣= 0.8. An inspection of Fig.3, where the time variations of the PDF W共s,t兲 moments are compared, demonstrates again that data from the numeric experiment and the analytic theory agree well within the lim-its of statistical accuracy共which is worse for N=256兲. Not surprisingly, the time scale of the MC results does not coin-cide with that of Eqs.共5兲 and 共6兲 since in the latter we have set K␣, which fixes the time scale, equal to unity. Closer examination of Fig.3shows that the resetting of the gener-alized diffusion coefficient as K␣⯝共80兲0.8⯝33.3 enables to superimpose the results of theoretical calculation and MC data.
In summary, we have shown unambiguously that the translocation dynamics of a polymer chain threading through a nanopore is anomalous in its nature. We have succeeded in calculating the anomalous exponent␣= 2 /共2+ 2 −␥1兲 from simple scaling arguments, and embedded it in the fractional diffusion formalism. We derived exact analytic expressions for the translocation time probability distribution as well as for the moments of the translocation coordinates, which are shown to agree well with our MC simulation data. The present treatment can be readily generalized to account for a drag force on the chain, and results for this case will be reported in a separate presentation.
The authors are indebted to Burkhard Dünweg for stimu-lating discussion during this study. The permission to use the special package for computation of Mittag-Leffler functions
withMATHEMATICA, provided by Y. Luchko, is gratefully
ac-knowledged. The authors are indebted to SFB DFG625 for financial support during this investigation.
关1兴 A. Meller, J. Phys.: Condens. Matter 15, R581 共2003兲. 关2兴 W. Sung and P. J. Park, Phys. Rev. Lett. 77, 783 共1996兲. 关3兴 P. J. Park and W. Sung, J. Chem. Phys. 108, 3013 共1998兲. 关4兴 M. Muthukumar, J. Chem. Phys. 111, 10371 共1999兲. 关5兴 H. Risken, The Fokker-Planck Equation 共Springer-Verlag,
Ber-lin, 1989兲.
关6兴 J. Chuang et al., Phys. Rev. E 65, 011802 共2001兲.
关7兴 Y. Kantor and M. Kardar, Phys. Rev. E 69, 021806 共2004兲. 关8兴 M. Doi and S. F. Edwards, The Theory of Polymer Dynamics
共Clarendon, Oxford, 1986兲.
关9兴 R. Metzler and J. Klafter, Phys. Rep. 339, 1 共2000兲. 关10兴 K. Luo et al., J. Chem. Phys. 124, 034714 共2006兲. 关11兴 D. Panja et al., e-print arXiv:cond-mat/0610671. 关12兴 R. Metzler and J. Klafter, Biophys. J. 85, 2776 共2003兲. 关13兴 P. G. de Gennes, Macromolecules 13, 1069 共1980兲.
关14兴 P. G. de Gennes, Adv. Colloid Interface Sci. 27, 189 共1987兲. 关15兴 C. Vanderzande, Lattice Models of Polymers 共Cambridge
Uni-versity Press, Cambridge, U.K., 1998兲.
关16兴 H. W. Diehla and M. Shpot, Nucl. Phys. B 528, 595 共1998兲; R. Hegger and P. Grassberger, J. Phys. A 27, 4069共1994兲. 关17兴 M. N. Barber et al., J. Phys. A 11, 1833 共1978兲.
关18兴 P. Debnath et al., J. Chem. Phys. 123, 204903 共2005兲; W. Min
et al., Phys. Rev. Lett. 94, 198302共2005兲; S. C. Kou and X. S.
Xie, ibid. 93, 180603共2004兲.
关19兴 W. T. Coffey et al., The Langevin Equation 共World Scientific, London, 2004兲.
关20兴 R. Metzler and J. Klafter, Physica A 278, 107 共2000兲. 关21兴 R. C. Lua and A. Y. Grosberg, Phys. Rev. E 72, 061918
共2006兲.
关22兴 S. B. Yuste and K. Lindenberg, Phys. Rev. E 69, 033101 共2004兲.
关23兴 It comes as no surprize that prefactors less than unity appear in 具s典共⬁兲 and 具s2典共⬁兲. For example, 具s2典共⬁兲 =兰0Ns2W共s,⬁兲ds/兰
0
NW共s,⬁兲ds=2N2, where 0⬍⬍1, and we have used the well-known mean value theorem for integration 共see any integral calculus textbook兲. This also does not contra-dict the condition that fixes the “experimental” average first-passage time : 具s2典共t=兲⬃N2, i.e., ⬃N2/␣ 关21兴. From the physical perspective, this means that the particular observation where s2共t=兲=N2is a relatively rear event in the total sam-pling.
关24兴 A. Milchev, K. Binder, and A. Bhattacharya, J. Chem. Phys.
121, 6042共2004兲.
关25兴 R. Gorenflo, J. Loutchko, and Yu. Luchko, Fractional Calculus Appl. Anal. 5, 491共2002兲.
FIG. 3.共Color online兲 Variation of the first and second moments of the PDF W共s,t兲 with time for chain lengths N
= 16, 32, 64, 128, 256. 共a兲 Log-log plot of the first moment 具s典 vs time t from a MC simulation共big symbols兲 and from Eq. 共5兲 共small
symbols兲. The dashed line denotes t␣/2with a slope of 0.4.共b兲 The same as in共a兲 but for 具s2典−具s典2. Analytical data are obtained from Eq.共6兲. The dashed line has a slope of 0.8.
DUBBELDAM et al. PHYSICAL REVIEW E 76, 010801共R兲 共2007兲
RAPID COMMUNICATIONS
