Pore-polymer interaction reveals non-universality in forced polymer translocation
We present a numerical study of forced polymer translocation by using two separate pore models. Both of them have been extensively used in previous forced translocation studies. We show that variations in the pore model affect the forced translocation characteristics significantly in the biologically relevant pore force, i.e. driving force, range. Details of the model are shown to change even the obtained scaling relations, which is a strong indication of strongly out-of-equilibrium dynamics in the computational studies which have not yet succeeded in addressing the characteristics of the forced translocation for biopolymers at realistic length scale.
Polymer translocation has been under intensive research for the past decade due to its relevance for e.g. ultra-fast DNA sequencing [1, 2, 3, 4, 5] and many biological processes . Most of the computational research on forced translocation, where the polymer is driven through a nanoscale pore by a potential, has dealt with fairly weak pore potential or force and assumed close-to-equilibrium dynamics. Recently, hydrodynamics has been shown to significantly affect forced translocation for experimentally and biologically relevant force magnitudes [7, 8]. The pore geometry has been noted to have effect on the experimental translocation process (see e.g. ). In addition, the effect of the pore model has been discussed in [9, 10]. However, the significance of the pore model in the forced polymer translocation has not yet been determined.
We have recently made a comparison of the unforced and forced translocation . In the first case, the process was seen to be close to equilibrium, as expected, and accordingly its dynamics turned out to be robust against variations in the computational model. However, in the second case the observed out-of-equilibrium characteristics for biologically and experimentally relevant pore force magnitudes, made the dynamics more sensitive to differences in the computational models [8, 12]. Accordingly, we expect the pore model to play a significant role in the forced translocation case .
So far, we have consistently used a cylindrical pore implemented by a potential symmetrical around the pore axis [8, 12, 11]. However, this differs from the typical pore implementation where the pore is formed by removing one particle from a wall consisting a monolayer of immobile, point-like particles [13, 14, 15, 16, 17].
A notable difference is that the latter pore implementation results in a driving potential that is inhomogeneous in the direction of the pore axis. In the case of lattice Boltzmann simulations these bead-pore models include both the square [16, 17], and cylindrically  shaped pores. Typically, the pore-polymer interaction is solely repulsive but the effect of attractive interaction has been investigated in a 2D model .
In order to investigate how susceptible translocation dynamics is to variations in the pore model we have implemented in our computational translocation model also the bead pore used in . We compare translocation processes for the bead pore and our original cylindrical pore. The model dependency is expected to be strongest in the pore region under strong forcing. We also try to evaluate the effect of the pore model for moderate forcing.
When the force bias inside the pore is small, the velocities involved are small enough for the hydrodynamics to be neglected [8, 11], and the Langevin dynamics (LD) is a valid choice for the computational model. Assuming that the polymer remains close to equilibrium throughout the translocation, the control parameter may be defined as the ratio of the total external force and the solvent friction, as done by Luo et al. in . Then the translocation dynamics would depend only on . They reported this to be the case even for a long polymer chain driven through the pore with a significant value of total pore force . The scaling of the average translocation time with respect to the chain length, would then be determined solely by . This is in clear contradiction with our earlier observation that the scaling exponent increases with the force due to crowding of the monomers on the trans side [8, 12]. In addition, the effect of the pore model along with other model-dependent factors would in this case be negligible .
The increase of with force cannot hold asymptotically as pointed out in , since ultimately very long polymers would translocate faster with a smaller pore force. From Fig. 1 b) in  it can be evaluated that this unphysical condition would be seen for for the pore force magnitudes and . This is comparable to actual DNA lengths of that are far beyond polymer lengths of used in simulations. The translocation dynamics then would have to be different for , which is in keeping with our earlier qualitative description of the forced translocation [8, 12]. The polymers were seen to be driven out of equilibrium throughout the translocation already for modest pore force magnitudes. The continually increasing drag force on the cis and crowding on the trans side were seen to change the force-balance condition, which is directly reflected on the translocation dynamics. This force balance is sure to be different for very long polymers.
The suggested control parameter , if applicable, would imply that the forced translocation could be completely characterized through varying it, which again is in stark contrast to our observation that the forced translocation is a non-universal, out-of-equilibrium process for biologically relevant magnitudes of pore force [8, 12]. Therefore, it is essential to determine if this contradiction arises from using different pore models.
This paper is organized such that we first describe the polymer model in Section II and the translocation models in Section III. The relation between computational and physical parameters is discussed in Section IV. The results are presented and discussed in Section V. We conclude the paper by summarizing our findings in Section VI.
Ii Polymer model
The standard bead-spring chain is used as a coarse-grained polymer model with the Langevin dynamics. In this model adjacent monomers are connected with an-harmonic springs, described by the finitely extensible nonlinear elastic (FENE) potential,
where is the length of an effective bond and is the maximum bond length. The Lennard-Jones (LJ) potential
is used between all beads of distance apart. The parameter values were chosen as , and . The used LJ potential mimics good solvent.
Iii Translocation models
To perform simulations of polymer translocation in 3D we use our translocation model based on Langevin dynamics . Hence, the time derivative of the momentum of bead reads as
where , , , and are the friction constant, the momentum, random force of the bead , and the external driving force, respectively. is constant but exerted only inside the pore. For unforced translocation . and are related by the fluctuation-dissipation theorem. The dynamics was implemented by using velocity Verlet algorithm .
In our model the wall containing the pore was implemented as a surface on which no-slip boundary conditions are applied to the polymer beads. The two different computational pore models to be compared can be described as follows (see Fig. 1):
Cylindrical pore. The pore, aligned with the -axis, is of diameter and length , where is the Kuhn length of the model polymer, and is either or . The pore is implemented by a cylindrically symmetric damped harmonic potential that pulls the beads toward the pore axis passing through the middle of the pore in the -direction. The pore is long, so the total pore force is taken as , where is the external driving force applied on each polymer bead inside the pore.
Bead pore. An octagonal composition of eight immobile particles was used as the bead pore model, as depicted in Fig. 1. This pore model includes a region where a polymer bead can reside without interacting with the pore beads, see Fig. 1. The effective pore length is approximately , so . Unlike in our model, in  also the wall was made of immobile particles. We verified that this has no effect on the translocation dynamics by reproducing essentially identical results on for the force magnitudes that were reported in .
In our previous Langevin dynamics simulations with the cylindrical pore [8, 11], parameter values , , and were used for the friction constant, the polymer bead mass, and the temperature in reduced units, respectively. Hence, for long times the one-particle self-diffusion constant was obtained from Einstein’s relation as . Time steps of and were (previously) used in the forced and unforced simulations, respectively. In this paper, to compare the two pore models, we have used , and as in , unless otherwise noted. Here, the time step is typically .
The number of beads is odd for polymers initially placed halfway inside the pore and even for polymers having initially only the first bead(s) inside the pore, see Fig. 1. Before allowing translocation a polymer is let to relax for longer than its Rouse relaxation time. We register events when a segment in the pore is replaced by the segment or . The polymer is considered translocated, when it has completely exited the pore to the trans side. Exit to either trans or cis side is regarded as an escape from the pore.
Iv Relating the computational and the physical pore force
The external driving force, called the total pore force, , where is the number of polymer beads inside the pore and the -directional driving force excerted on each bead inside the pore, will be taken as the effective pore force in our model. The thermal energy is set to in reduced units and the length scale is set by the polymer bond length . The magnitude of the effective pore force is then determined with respect to thermal fluctuations by comparing with . It is noteworthy, however, that since the pore length is of the order of the polymer bond length, there exists an inevitable computational artefact due to the discrete polymer model. The external force is excerted on beads residing in the pore, whose number does not stay constant during translocation. In addition this number is different for different pore lengths, so the pore force does not increase strictly linearly with the pore length. These effects are, however, mitigated by the fact that averaged over translocation time is constant.
Correspondence between computational and physical length scale can be established by taking the polymer bond length as the Kuhn length for the physical polymer. In SI units the bond length for our FJC model polymer can be obtained as , where is the persistence length, for a single-stranded (ss) DNA . The total pore force in SI units, , is then obtained from the dimensionless total pore force, , by relating . Here, is the Boltzmann constant, and the physical temperature is taken to be . Thus the effective pore force of corresponds to for a ssDNA. To relate to experiments, in the -HL pore a typical pore potential of would correspond to when Manning condensation leading to drastic charge reduction is taken into account [22, 23]. The used computational pore force magnitudes are thus in the physically relevant range .
v.1 Forward transfer probability
Previously , we have shown that during the unforced translocation the polymer remains close to equilibrium so that the forward transition probabilities
derived from the free energy apply in 3D, where . Including the pore force in the calculation only introduces a prefactor to the above not changing its form, i.e. dependence on . For this to hold, the translocating polymer has to remain close to equilbrium. Indeed, except for a small increase due to the force term, the shapes of the measured for the unforced and forced translocation with the total pore force are seen to be similar, cf. Fig. 2. for the two pore models do not differ appreciably for this small pore force. Hence we conclude that the transition probabilities can be obtained from equilibrium framework for pore force magnitudes of the order of .
For the bead pore already for the force , the curve differs from the equilibrium shape, as shown in Fig. 2. Hence, the effect of the pore force can no more be taken as a small perturbation to equilibrium dynamics. In keeping with this, the transition probabilities for are seen to depend on the polymer’s initial position, here either on the cis side or halfway through the pore. In other words, the transition probability from the current state to the next depends on the path through which this state was arrived at. The polymer cannot then have relaxed to equilibrium between previous transitions.
Using the cylindrical pore, the probability is seen to be almost constant and close to one for all with , see Fig. 2. This is in accordance with our previous results that the polymer was seen to almost always translocate with , which means that also in this case would be close to one . However, for the bead pore deviates significantly from for the cylinder pore when , see Fig. 2. Thus it is evident that the pore force magnitudes in the two model pores do not have an exact correspondence for large pore force magnitudes.
The pore force at which translocation is seen to be a strongly out-of-equilibrium process is surprisingly low for both pore models. It may be compared with the random force in the one-particle Langevin equation satisfying fluctuation-dissipation theorem in 3D, for which . With the parameter values and used in our simulations, this yields .
Already at fairly moderate pore force values, a local force balance governs the translocation dynamics. This was first proposed by Sakaue , and demonstrated in simulations by us [8, 12]. Considering a concept of ‘mobile beads’ near the pore, we obtained and described qualitatively the scaling relation , where the parameter taking into account the varying number of mobile beads diminishes toward zero with increasing pore force [8, 12]. Sakaue has recently given a fairly quantitative out-of-equilibrium formulation for the scaling . However, neither description takes into account the crowding of monomers on the trans side shown to be significant for these pore force magnitudes .
Here, we present results from simulations of the polymer translocation covering the relevant range of the pore force. When , attains the equilibrium shape and when , is close to one.
v.2 The proposed control parameter
Using the bead pore, we first measured the scaling exponent for different values of pore force while keeping the proposed control parameter fixed. For the parameter pair values , , and were used, see Fig. 3 (a). Similar average translocation times were obtained for large pore force magnitudes, . However, differs appreciably for , resulting in a different value for the scaling exponent .
The measured average waiting (or state transition) times are compared in Fig. 3 (b). The shown three plots for are for the above-mentioned values of . In the limit , the average state transition times were seen to be similar. This is in accord with the transition times in the unforced translocation obeying a close-to-Poissonian distribution . When the friction of the solvent dominates over the stochastic term in the Langevin equation, Eq. (3). The transition time profiles for the bead pore are qualitatively similar with those obtained for the cylindrical pore with different parameters .
From the different scaling of the translocation times obtained for large and small pore force magnitudes and the different average state transition times we conclude that which was kept fixed cannot be a universal control parameter for forced translocation regardless of the pore model. This is in keeping with our previous finding that the forced translocation is a highly out-of-equilibriun, non-universal process for the biologically and experimentally relevant pore force range [8, 12].
v.3 Pore model matters in forced translocation
Originally, the need to compare the cylindrical and bead pore models arose from the differences in the scaling of the translocation time with the polymer length in forced translocation obtained using these two pore models. The translocation times averaged over at least and at most runs are shown in Fig. 4 (a) as functions of . The total pore force is constant, , and chosen to be within the experimentally relevant range . For the friction parameter three values, , , and , are used. For the two separate pore models we find clear differences in both the absolute value of and its apparent scaling. Regardless of the value for , we obtain a smaller scaling exponent for the cylindrical than for the bead pore. In addition, increasing from to increases for both pore models, see inset in Fig. 4 (b). This in part explains the different results obtained with the two pores as the effective friction experienced by a chain is different inside them. This difference is naturally enhanced at higher velocities.
Our previous results were obtained with the polymer bead mass due to the requirement of the stochastic rotation dynamics that the polymer beads should be heavier than solvent beads [8, 12]. Therefore, we checked how the bead mass affects the characteristics of the translocation time. Increasing the bead mass from to is seen not only to increase the average translocation time but also to change its apparent scaling with polymer length , see Fig. 4 (a). To be certain, we checked this using two independent Langevin algorithms.
In accordance with our previous findings for cylindrical pore, was seen to increase with pore force also for the bead pore. For , the scaling exponent (Fig. 4 (b)), (Fig. 3 (a)), (Fig. 4 (a)), and (not shown) were obtained, respectively. Surprisingly, for the cylindrical pore , which increased substantially with increasing when [8, 12], did not show such strong tendency when . Changing seems then to change even the qualitative characteristics of forced polymer translocation.
The change of in the scaling due to a change in is hardly surprising. However, changing with is a bit more subtle. It suggests that we are actually looking at the forced translocation in a continually changing or transient stage. This is in accord with our previous finding that for polymers of lengths that are used in simulations, , the number of moving polymer beads changes throughout the translocation and thus continually alters the force balance condition resulting from the forced translocation taking place out of equilibrium . Also the crowding of the polymer beads on the trans side, whose relaxation toward equilibrium is slower than the translocation rate changes the force balance continually. Hence, the simulated forced translocation processes do not reveal the asymptotic () characteristics for experimentally relevant pore force magnitudes. Then reporting scaling exponents for such forced translocation seems unwarranted. The remaining question of interest then is, why does the forced translocation exhibit the scaling albeit with varying exponents?
v.4 Translocation at low pore force
For simulations at low pore force the chains were initially placed halfway through the pore unlike at large pore force where the polymer was placed so that only its end was inside the pore. At the very low force of the pore model was found to affect only slightly the absolute value of , see Fig. 4 (b). For the bead pore model chains of length had a finite probability to escape to the cis side. Unlike for the bead pore, even some of the beads long polymers escaped to the cis side for the cylindrical pore. This is probably due to the cylindrical pore aligning the polymer chain toward the pore axis in the middle thus effectively reducing the friction between the polymer and the pore.
For the pore force we obtain , see Fig. 4 (b), which is the expected exponent for unforced translocation, , since has been measured for the swelling exponent in our model . Approaching the unforced case makes the translocation dynamics increasingly robust to variations in the model, which is a natural consequence of the small pore force not dominating over entropic forces. We regard the pore force magnitude already small enough for the translocating polymer to remain close to equilibrium, see Fig. 2. However, we did not obtain a clear scaling for for either pore models, but close to the unforced translocation value was obtained for chains shorter than while a lower value of was obtained for longer chains. This applies for both and , see Fig. 4 (b). We expect the lower value of to be (see ), although a wider range of would be needed to confirm this. This apparent cross-over behavior for low pore force is outside the scope of the present paper and calls for a separate, more thorough investigation.
For the presumably small pore force we obtained . This differs considerably from reported in . From Fig. 2 it can be seen that for clearly deviates from the close-to-unforced form of for . Together with the low value of this observation suggests that is large enough to drive the polymer more and more out of equilibrium as the translocation proceeds. The different values for here and in  also support this view, since the two simulations started from different initial positions. When the relaxation time of the polymer toward thermal equilibrium is slower than the translocation time, the process possesses memory, or correlation over time, i.e. the memory function . Thus, a polymer starting half-way through the pore is in a different state than a polymer that has arrived at this position and started from another initial position.
In summary, we have simulated forced polymer translocation in 3D by using Langevin dynamics. We have implemented two pore models in our algorithm, (i) a pore surrounded by eight immobile beads, which we call bead pore and (ii) a cylindrical pore where a damped harmonic potential confines the beads inside the pore region. The present study compares the effect these two pores have on the forced translocation.
We measured the forward transition probabilities to determine the range of pore force that would be sufficiently large to include cases where polymer translocation takes place close to and strongly out of equilibrium. We found that the polymer remains close to equilibrium when the total pore force . Here the translocation dynamics was found to be robust to variations in the translocation model, such as the details of the pore. was found to deviate significantly from the close-to-equilibrium form already for a pore force as small as . The polymer was found to be driven far from equilibrium when .
For small pore force magnitudes the forced translocation processes are identical for the two pore models. However, the translocation characteristics were found to be increasingly model dependent when the pore force is increased. This is a natural consequence of the dynamics of the forced translocation being determined by a continually changing force-balance condition when the pore force is large enough, i.e. (). Accordingly, it seems that universal exponents for the forced translocation cannot be found in the biologically relevant pore force regime. In addition, attempts to define a control parameter whose magnitude would consistently determine these exponents would seem futile, which we showed for one proposed candidate by using the bead pore. Qualitatively the forced translocation exhibited similar characteristics with both pore models, most notably the increase with the pore force of the exponent determining the relation between the average translocation time and the polymer length. In our view this is another indication of the continually changing force-balance condition governing the highly non-equilibrium forced translocation process.
Acknowledgements.One of the authors (V.V.L.) thanks Dr. K. Luo for email exchanges. This work has been supported by the Academy of Finland (Project No. 127766). The computational resources of CSC-IT Centre for Science, Finland, are acknowledged.
- Kasianowicz et al.  J. J. Kasianowicz, E. Brandin, D. Branton, and D. W. Deamer, Proc. Natl. Acad. Sci. U.S.A. 93, 13770 (1996).
- Storm et al  W. Storm et al, Nano Lett. 5, 1193 (2005).
- Li et al.  J. Li, M. Gershow, D. Stein, E. Brandin, and J. A. Golovchenko, Nature Materials 2, 611 (2003).
- Meller et al.  A. Meller, L. Nivon, and D. Branton, Phys. Rev. Lett. 86, 3435 (2001).
- Zwolak and Di Ventra  M. Zwolak and M. Di Ventra, Rev. Mod. Phys. 80, 141 (2008).
- Alberts et al.  B. Alberts et al., Molecular Biology of the Cell (Garland Publishing, New York, 1994).
- Fyta et al. [2008a] M. Fyta, S. Melchionna, S. Succi, and E. Kaxiras, Phys. Rev. E 78, 036704 (2008a).
- Lehtola et al.  V. V. Lehtola, R. P. Linna, and K. Kaski, EPL 78, 58006 (2009).
- Lubensky and Nelson  D. K. Lubensky and D. R. Nelson, Biophys. J. 77, 1924 (1999).
- Slonkina and Kolomeisky  E. Slonkina and A. B. Kolomeisky, J. of Chem. Phys. 118, 7112 (2003).
- Lehtola et al.  V. V. Lehtola, R. P. Linna, and K. Kaski, Phys. Rev. E 81, 031803 (2010).
- Lehtola et al.  V. V. Lehtola, R. P. Linna, and K. Kaski, Phys. Rev. E 78, 061803 (2008).
- Luo et al.  K. Luo, S. T. T. Ollila, I. Huopaniemi, T. Ala-Nissila, P. Pomorski, M. Karttunen, S.-C. Ying, and A. Bhattacharya, Phys. Rev. E 78, 050901(R) (2008).
- Luo et al.  K. Luo, T. Ala-Nissila, S.-C. Ying, and R. Metzler, EPL 88, 68006 (2009).
- Gauthier and Slater  M. G. Gauthier and G. W. Slater, Eur. Phys. J. E 25, 17 (2008).
- Fyta et al. [2008b] M. Fyta, E. Kaxiras, S. Melchionna, and S. Succi, Comp. Sci. & Eng. 2, 20 (2008b).
- Fyta et al.  M. Fyta, S. Melchionna, E. Kaxiras, and S. Succi, Multiscale Modeling and Simulation 5, 1156 (2006).
- Bernaschi et al.  M. Bernaschi, S. Melchionna, S. Succi, M. Fyta, and E. Kaxiras, Nanoletters 8, 1115 (2008).
- Luo et al.  K. Luo, T. Ala-Nissila, S.-C. Ying, and A. Bhattacharya, Phys. Rev. Lett. 99, 148102 (2007).
- van Gunsteren and Berendsen  W. van Gunsteren and H. Berendsen, Mol. Phys. 34, 1311 (1977).
- Tinland et al.  B. Tinland, A. Pluen, J. Sturm, and G. Weill, Macromol. 30, 5763 (1997).
- Meller  A. Meller, J. Phys. Condens. Matter 15, R581 (2003).
- Sauer-Budge et al.  A. F. Sauer-Budge, J. A. Nyamwanda, D. K. Lubensky, and D. Branton, Phys. Rev. Lett. 90, 238101 (2003).
- Sakaue  T. Sakaue, Phys. Rev. E 76, 021803 (2007).
- Sakaue  T. Sakaue, Phys. Rev. E 81, 041808 (2010).
- Kantor and Kardar  Y. Kantor and M. Kardar, Phys. Rev. E 69, 021806 (2004).