Multiscale Coupling of Molecular Dynamics and Hydrodynamics: Application to DNA Translocation through a Nanopore

. We present a multiscale approach to the modeling of polymer dynamics in the presence of a ﬂuid solvent. The approach combines Langevin Molecular Dynamics (MD) techniques with a mesoscopic Lattice-Boltzmann (LB) method for the solvent dynamics. A unique feature of the present approach is that hydrodynamic interactions between the solute macromolecule and the aqueous solvent are handled explicitly, and yet in a computationally tractable way due to the dual particle-ﬁeld nature of the LB solver. The suitability of the present LB-MD multiscale approach is demonstrated for the problem of polymer fast translocation through a nanopore. We also provide an interpretation of our results in the context of DNA translocation through a nanopore, a problem that has attracted much theoretical and experimental attention recently.

are changing the scene. Minimal forms of the Boltzmann equation can be designed on the lattice, which quantitatively describe the behavior of fluid flows in a way that is often computationally more advantageous than the continuum approach based on the Navier-Stokes equations. Moreover, lattice kinetic theory has proven capable of dealing with complex flows, such as flows with phase transitions and strong heterogeneities, for which continuum equations are exceedingly difficult to solve, if at all known (for a recent review see [6]). These advances have opened the road to developing new mesoscopic multiscale solvers [7]. The present work provides a successful implementation of such an approach. We will focus on the coupling of a mesoscopic fluid solver, the lattice Boltzmann method, with simulations at the atomistic scale employing explicit molecular dynamics. A unique feature of our approach is the dual nature of the mesoscopic kinetic solver, which propagates coarse-grained information (the single-particle Boltzmann probability distribution), along straight particle trajectories. This dual field/particle nature greatly facilitates the coupling between the mesoscopic fluid and the atomistic levels, both on conceptual and computational grounds.
The paper is organized as follows. In Section §2 we present the basic elements of the multiscale methodology, namely the Lattice Boltzmann treatment of the fluid solvent, and its coupling to a Molecular Dynamics simulation of the solute biopolymer. In Section §3, we present an application of this multiscale methodology to the problem of long polymer translocation through a nanopore; in particular, we analyze in detail the role of hydrodynamics in accelerating the translocation process. In Section §4 we elaborate on the relevance of our results to the problem of DNA translocation, which has attracted much theoretical and experimental attention recently. We conclude in Section §5 with general remarks and outlook for future extensions.

Lattice-Boltzmann -Molecular-Dynamics multiscale methodology.
We consider the generic problem of tracing the dynamic evolution of a polymer molecule interacting with a fluid solvent. This involves the simultaneous interaction of several physical mechanisms, often acting on widely separate temporal and spatial scales. Essentially, these interactions can be classified in three distinct categories as solute-solute, solvent-solvent and solvent-solute. The first category includes the conservative many-body interactions among the single monomers in the polymer chain. Being atomistic in nature, these interactions usually set the shortest scale in the overall multiscale process. They are typically handled by Molecular Dynamics techniques for constrained molecules. The second category, the solvent-solvent interactions, refer to the dynamics of the solvent molecules, which are usually dealt with by a continuum fluid-mechanics approach; in the present work these will be described by the mesoscopic Lattice Boltzmann equation. The second and third category have also been handled by simulating the solvent explicitly via molecular dynamics, implicit solvent particles via Brownian dynamics including hydrodynamic interactions, or solving the corresponding Fokker-Planck equation [8]. Finally, the solvent-solute dynamics will be treated by augmenting the molecular dynamics side with dissipative fluid-molecule interactions (Langevin picture) and including the corresponding reaction terms in the fluid-kinetic equations.
2.1. Atomistic dynamics. We consider a polymer consisting of N monomer units (also referred to as beads). The polymer is advanced in time according to the following set of Molecular Dynamics-Langevin equations for the bead positions r p and velocities v p : where we distinguish four types of forces: The first term represents the conservative bead-bead interactions through a potential which we will take to have the standard 6 − 12 Lennard-Jones form, truncated at a distance of r=2 1/6 σ [9]. This was combined with a harmonic part to account for the energy cost of distorting the angular degrees of freedom, with φ the relative angle between two consecutive bonds. Torsional motions will not be included in the present model, but can easily be incorporated if needed.
We consider next the solute-solvent interactions. The second term on the righthand-side of Eq.(2.1) represents the mechanical friction between the single bead and the surrounding fluid, v p being the bead velocity and u p the fluid velocity evaluated at the bead position. In addition to mechanical drag, the polymer feels the effects of stochastic fluctuations of the fluid environment, through the random term, ξ p , a Gaussian noise obeying the fluctuation-dissipation relations: where V is the volume of the cell to which beads p and q belong. Finally, λ p ∂ rp κ p is the reaction force resulting from N − 1 holonomic constraints for molecules modelled with rigid covalent bonds: r 0 being the prescribed bond length, and {λ p } is the set of N − 1 Lagrange multipliers conjugated to each constraint. The usage of constraints instead of flexible bond lengths makes it possible to eliminate unimportant high-frequency intra-molecular motion which would render the underlying LB propagation prone to numerical instabilities. In this way, the time-step of the Molecular Dynamics part can be increased by about one order of magnitude, as much as the overall efficiency of the LBMD method, as we shall discuss in section §2.5. Finally, in order to avoid spurious dissipation, the bead velocities are required to be strictly orthogonal to the relative displacements. Given the second order atomistic dynamics, the velocities must obey the independent constraints: The constraints (2.9), (2.10) are enforced over positions and momenta separately via the SHAKE [10] and the RATTLE algorithms [11]. The implementation of these constraints requires the iterative solution of the system of equations (2.9)-(2.10), typically accomplished via standard Newton-Raphson techniques.

2.2.
Fluctuating Lattice Boltzmann method. The Lattice Boltzmann equation is a minimal form of the Boltzmann kinetic equation in which all details of molecular motion are removed except those that are strictly needed to recover hydrodynamic behavior at the macroscopic scale (mass-momentum and energy conservation). The result is an elegant equation for the discrete distribution function f i ( x, t) describing the probability to find a particle at lattice site x at time t with speed v = c i . More specifically, since we are dealing with nanoscopic flows, in this work we shall consider the fluctuating Lattice Boltzmann equation which takes the following form: where f i ( x, t) represents the probability of finding a fluid particle at spatial location x and time t with discrete speed c i . The particles can only move along the links of a regular lattice defined by the discrete speeds, so that the synchronous particle displacements ∆ x i = c i ∆t never take the fluid particles away from the lattice. For the present study, the standard three-dimensional 19-speed lattice is used [4] (see Figure 1). The right hand side represents the effect of intermolecular solvent-solvent collisions, through a relaxation toward local equilibrium, f eq i , typically a second order (low-Mach) expansion in the fluid velocity of a local Maxwellian with speed u: where β = m f /k B T f is the inverse fluid temperature (with k B the Boltzmann constant), w i a set of weights normalized to unity, and I is the unit tensor in configuration space. The relaxation frequency ω controls the fluid kinematic viscosity ν through the relation: where c s is the sound speed in the solvent [7]. Knowledge of the discrete distributions f i allows the calculation of the local density ρ, flow speed ρ u and momentum-flux tensor ↔ P , by a direct summation upon all discrete distributions: The diagonal component of the momentum-flux tensor gives the fluid pressure, while the off-diagonal terms give the shear-stress. Unlike in hydrodynamics, both quantities are available locally and at any point in the simulation.
Thermal fluctuations are included through the source term F i which reads as follows (index notation) where F (2) is the fluctuating stress tensor (a 3 × 3 stochastic matrix). Consistency with the fluctuation-dissipation theorem at all scales requires the following conditions where ∆ abcd is the fourth-order Kronecker symbol [12]. F (3) is related to the fluctuating heat flux and g iabc is the corresponding basis in kinetic space, essentially a third-order Hermite polynomial (full details are given in [13]).
The polymer-fluid back reaction is described through the source term S i , which represents the momentum input per unit time due to the reaction of the polymer on the fluid population f i : where D(x) denotes the mesh cell to which the p th bead belongs. The quantities on the left hand side in the above expression have to reside on the lattice nodes, which means that the frictional and random forces need to be extrapolated from the particle to the grid location.
The use of a LB solver for the fluid solvent is particularly well suited to this problem because of the following reasons: i) Free-streaming proceeds along straight trajectories. This is in stark contrast with hydrodynamics, in which fluid momentum is transported by its own space-time varying velocity field. Besides securing exact conservation of mass and momentum of the numerical scheme, this also greatly facilitates the imposition of geometrically complex boundary conditions.
ii) The pressure field is available locally, with no need of solving any (computationally expensive) Poisson problem for the pressure, like in standard hydrodynamics. iii) Unlike hydrodynamics, diffusivity is not represented by a second-order differential operator, but it emerges instead from the first-order LB relaxation-propagation dynamics. The result is that the kinetic scheme can march in time-steps which scale linearly, rather than quadratically, with the mesh resolution. This facilitates highresolution down-coupling to atomistic scales. iv) Solute-solute interactions preserve their local nature since they are explicitly mediated by the solvent molecules through direct solvent-solute interactions. As a result, the computational cost of hydrodynamic interactions scales only linearly with the length of the polymer (no long-range interactions). v) Since all interactions are local, the LB scheme is ideally suited to parallel computing.
It is worth mentioning that more advanced Lattice Boltzmann models [14,15] could equally well be coupled to the atomistic dynamics. The numerical solution of the stochastic equations is performed by means of a modified version of the Langevin Impulse propagation scheme, derived from the assumption that the systematic forces are constant between consecutive time steps [16]. The propagation of the unconstrained dynamics proceeds according to the scheme [17]r where C(dt) is an array of 3N gaussian random variables with zero mean and variance kB T m (e 2γdt − 1), and {r} represent temporary positions. The propagator (2.20) is particularly suitable for our purposes since it is second order accurate in time and robust, that is, it reduces to the symplectic Verlet algorithm for γ → 0. Moreover, at variance with the original Langevin Impulse scheme, the modified propagator allows for an unambiguous definition of velocities, which are needed to couple the polymer to the hydrodynamic field of the surrounding solvent. The particle positions and velocities corrected via the SHAKE and RATTLE algorithms read For consistency, in considering the momentum exchange with the solvent the corrected velocities appear in the friction forces. The MD cycle is repeated M times, with the hydrodynamic field frozen at time t n = n∆t.
2.4. Spatial exchange. The transfer of spatial information from/to grid to/from particle locations is performed at each LB time-stamp t n = n∆t. To this purpose,on account of its simplicity, a simple nearest grid point (NGP) interpolation scheme is Transfer of spatial information (a) from grid to particle, and (b) from particle to grid. Black spheres denote beads, while in white are the lattice sites.
used (see Fig.2). Momentum conservation was checked to hold up to six digits. With reference to a time slice t n = n∆t, the pseudo-algorithm performing a single LB time-step, reads as follows 1. Interpolation of the velocity: u( x) → u p 2. For m = 1, M : Advance the molecular state from t to t + dt 3. Extrapolation of the forces: F p → F ( x) 4. Advance the Boltzmann populations from t to t + ∆t This time-marching can be formally represented by an operator-splitting multistep time procedure for two coupled kinetic equations describing the dynamic evolution for the fluid and the polymer distribution functions, respectively [18]. It is worth emphasizing that, while LB and MD with Langevin dynamics have been coupled before, notably in the study of single-polymer dynamics [19], to the best of our knowledge, this is the first time such that coupling is put in place for long molecules of biological interest.

Efficiency considerations.
The total cost of the computation scales roughly like where t LB is the CPU time required to update a single LB site per timestep and t MD is the CPU time to update a single bead per timestep, V is the volume of the computational domain in lattice units and N is the number of polymer beads, with M the LB-MD time-step ratio. Finally, N LB is the number of LB timesteps. In the above equation, t MD includes the overhead of LB-MD coupling. Note that t MD is largely independent of N because i) the LB-MD coupling is local, ii) the forces are short ranged and iii) the SHAKE/RATTLE algorithms are empirically known to scale linearly with the number of constraints.
Regarding the cost of the LB section, this is known to scale linearly with the volume occupied by the solvent. For the case where polymer concentration is kept constant, the volume needed to accommodate a polymer of N beads should scale approximately as N 1.8 ; however, for translocation studies such as those discussed later in this paper, we shall consider a box of given volume, independently on the polymer length.
From the above expression it is clear that M should be chosen as small as possible, consistent with the requirement of providing a realistic description of the polymer dynamics. In the present simulation we typically choose M between 5 and 20, depending on the parameters of the simulation, particularly the temperature. This means that we are taking the LB representation close to the molecular scale. We will return to this important issue in the quantitative discussion of the physical application.
A tentative estimate of the computational cost proceeds as follows: Assuming 250 flops/site/LB-step and 2500 flops/bead/MD-step (including the LB-MD coupling overhead), and an effective processing speed of 100 Mflop/s, the evolution over 30, 000 LB steps=150, 000 MD steps of a typical 80×40×40 grid and 400 beads set-up, would take about: which is in reasonable agreement with the simulation time observed with the present version of the code (7hrs), including the relative MD/LB cost (≃ 1 : 4).
We wish to emphasize that the key feature of the LB-MD approach, namely linear scaling of the CPU cost with the number of beads (at constant volume) is indeed observed. In fact, the execution times for 50, 100 and 400 beads are 0.433, 0.489, and 0.882 sec/step, respectively on a 2GHz AMD Opteron processor. By excluding hydrodynamics, these numbers become 0.039, 0.075, and 0.318 sec/step. It is worth mentioning that thus far, no effort has been directed to code optimization; it is quite possible that careful optimization may lower the execution time by an order of magnitude.
2.6. Validation tests. The static and dynamic behavior of the DNA chain obtained by our methodology has been compared to the scaling predictions for a single chain at infinite dilution. Given the structure factor S f (k) = 1 N i,j e ik·(ri−rj ) , standard theory predicts the scaling law S f (k) = N g(kR g ), where g(y) is a universal function and R g = (1/2N 2 ) i,j (r i − r j ) 2 is the gyration radius. For large k, the structure factor is independent of N , and it follows where experiments, theory and simulations agree on the scaling exponent value µ ≃ 0.584. The static scaling law is not affected by the presence of hydrodynamics. However, verification of the scaling law and attainment of the scaling regime for large enough chains is a good check for the correctness of our simulation scheme and for the subsequent validation of the hydrodynamic behavior.
The dynamic behavior of the chain is deeply affected by the presence of hydrodynamic interactions. The standard picture of polymer dynamics is based on the Rouse (no hydrodynamics) or Zimm (hydrodynamic) description in terms of an underlying gaussian chain. In this case, the chain intermediate scattering function whereg(y) is another universal function and D f is the center of mass diffusion constant. Having introduced the dynamic scaling exponent µ ′ via D f ∼ N µ ′ , the exponent x is found to be x = 2 + µ ′ /µ. According to Zimm theory, µ ′ = µ and x = 3 [19] while, according to Molecular Dynamics simulations, it appears that the actual value is somehow lower, i.e. x ≃ 2.9 [20].
We have considered a chain made of 30 monomers with bead-bead Lennard-Jones parameters taken from previous studies of chains simulated via Brownian dynamics [21] (σ = 0.65, ǫ/k B T = 1.0, bond length r 0 = 0.945) in a simulation box of edge 60. This choice was motivated to verify the range of scaling behavior as compared to previous numerical results. We have computed the structure factor, as reported in Fig.3(a), and observed that the scaling regime is clearly visible for 1 < k < 3 with an exponent equal µ = 0.58 ± 0.1. In the range of k vectors where the static scaling holds, the dynamic scaling has been checked by considering, for given values of the computed scattering function, the loci of points t(k, I S ) =g −1 (I S )k −x and by fitting via a power law curve (see ref. [20] for details). As illustrated in Fig.3(b), the resulting scaling exponent is found to be x = 2.9±0.1, in excellent agreement with the expected value and similar to previous simulation results on single polymers surrounded by a Lattice Boltzmann fluid [21]. Moreover, by applying a heuristic argument [20], we have verified that within the scaling regime, the finite size of the periodic box was not biasing the data.
3. Application: polymer translocation through nanopore. The scheme described above is general and applicable to any situation where a long polymer is moving in a solvent. This motion is of great interest for a fundamental understanding of polymer dynamics in the presence of the solvent. For example, the translocation of a polymer through a pore of very small size (of order the separation between monomers), is a process in which the coupling of the molecular motion to the solvent dynamics may be of crucial significance. In this section, we will therefore provide a detailed discussion of the polymer dynamics in the presence of a solvent for the example of translocation through a nanopore but without reference to a specific physical system. In the next section we explore the relevance of these results to DNA translocation through a nanopore.
3.1. Initial and Boundary conditions. The polymer is initialized via a standard self-avoiding random walk algorithm and further relaxed to equilibrium by standard Molecular Dynamics. The solvent is initialized with the equilibrium distribution corresponding to a constant density ρ 0 and zero macroscopic speed u = 0.
Boundary conditions for the fluid are periodic at inlet/outlet sections, and zerospeed at rigid walls, using the standard bounce-back rule [4]. For the polymer, periodicity is again imposed at inlet/outlet, whereas the interaction with rigid walls is han-dled by a Lennard-Jones potential with specific wall-polymer parameters σ wall = 1.5 and ǫ wall = 10 −3 in LB units. The connection between slip-flow at the wall and intermolecular solid-fluid interactions shall be the objects of future research.

Numerical set-up.
We consider a three-dimensional box of size N x h × N y h × N z h lattice units, with h the spacing between lattice points. We will take N x = 2N y , N y = N z ; the separating wall is located in the mid-section of the x direction, at x = hN x /2 with N x = 80. At t = 0 the polymer resides entirely in the right chamber at x > hN x /2. At the center of the separating wall, a square hole of side d hole = 2h is opened, through which the polymer can translocate from one chamber to the other. Translocation is induced by a constant electric field which acts along the x direction, and is confined in a rectangular channel of size 2h × h × h along the streamline (x direction) and cross-flow (y, z directions). The spatial coarse-graining is such that the presence of the solvent as well as electrostatic forces acting due to charges on the polymer are neglected altogether as being of secondary importance compared to hydrodynamics.
Here and throughout we work in lattice Boltzmann units, in which length and time are measured in units of the lattice spacing h = ∆x and time-step ∆t, respectively. Mass is defined as m = m LB m sol . The dimensionless mass m LB used in the simulations is set to unity, which means that mass is measured in units of the solvent mass m sol . This choice is not restrictive since the present approach is used to model incompressible flows in which density is a parameter which can be rescaled by any arbitrary factor. However, it is of some interest to estimate the number of solvent molecules represented by a single LB computational molecule, since the inverse of this number conveys a measure of the importance of statistical fluctuations at the scale of the lattice spacing ∆x. Let S N be this number, which will be defined as ρ sol h 3 ρLB m sol , where ρ LB is the dimensionless density used in the LB simulations. In order for the Boltzmann probability distribution to make sense as a statistical observable, N LB >> 1. For typical values of ρ sol = 1 gr/cm 3 , m sol ∼ 20 amu, ρ LB = 1 (which correspond to water), and h in the range 1 − 10 2 nm, this yields S N ∼ 10 4 − 10 6 . This shows that the neglect of many-body fluctuations inherent to the single-particle Boltzmann representation is still justified even at the nanoscopic scale of the lattice spacing.
We will focus here on the fast translocation regime, in which the translocation time t X is much smaller than the Zimm time, t Z , i.e. the typical relaxation time of the polymer towards its native (minimum energy, maximum entropy) configuration. Under fast-translocation conditions, the many-body aspects of the polymer dynamics cannot be ignored because different beads along the chain do not move independently. As a result, simple one-dimensional Brownian models do not apply [22]. In addition to many-body solute-solute interactions, the present approach also takes full account of many-body solute-solvent hydrodynamic interactions. The conditions for fast-translocation regime can be appraised as follows. The translocation time is estimated by equating the driving force, F pull , to the drag force exerted by a solvent with dynamic viscosity η on a polymer with radius of gyration R, F drag ∼ 6πηR 2 tX . This yields t X ∼ 6πηR 2 F pull . Since the Zimm time is given by t Z ∼ 0.4ηR 3 kB T , the fast-translocation condition t X ≪ t Z becomes: Our reference simulation is performed with F pull /m = 0.02∆x/∆t 2 and k B T /m = 10 −4 ∆x 2 /∆t 2 , with m the mass of one bead (monomer) of the polymer. The polymer length is in the range 20 ≤ N ≤ 400 beads. It can be readily checked that by assuming R ∼ N 0.6 our set of parameters falls safely within the fast translocation regime. However, for k B T /m = 10 −3 , F pull R/k B T is of the order of 10 2 − 10 3 which is much closer to breaking the above condition.
The main parameters of the simulation are (in LB units) σ = 1.8 and ǫ = 10 −4 for the Lennard-Jones potential. The bond length among the beads is set at r 0 = 1.2. According to these values, the Lennard-Jones time-scale, τ LJ = σ/ 2ǫ/m, is of the order of ∼ 100∆t. Thus, by choosing M = 5 as a time-gap factor, we obtain dt ∼ τ LJ /500, which is adequate for the resolution of the polymer dynamics. The solvent is set at a density ρ = 1, with a kinematic viscosity ν = 0.1∆x 2 /∆t and a damping coefficient γ = 0.1/∆t. The flexional rigidity κ for the angular potential between beads will be 10 −4 /rad. In order to resolve the structure of the solvent accurately on the atomistic scale [23], we should use a higher resolution of at least 3-4 orders of magnitude. This means resolving the radial structure of the pore, a task that can only be undertaken by resorting to parallel computing. It is nonetheless hoped, and verified a posteriori, that this artificial magnification does not affect adversely the most significant dynamical and statistical properties of the translocation process, by which we mean that eventually, the time-scale of the simulated process may not be the same as in the physical process of interest, but the simulated dynamics is related to the physical dynamics by a simple rescaling of the time variable.

Translocation time.
The most immediate quantity of interest in the translocation process is the dependence of the translocation time on the polymer length. This is usually expressed by a scaling law of the form where t 0 is a reference time-scale, formally the translocation time of a single monomer, and α a scaling exponent measuring the degree of competition (α > 1) / cooperation (α < 1) of the various monomers in the chain.
We first turn to the derivation of the scaling behavior of the translocation process in the case where hydrodynamic interactions are included. In order to take into account the statistical nature of the phenomenon, simulations of a large number of translocation events (100 up to 1000) for each polymer length were carried out. The ensemble of simulations is generated by different realizations of the initial polymer configuration. The duration histograms were constructed by cumulating all events for various lengths. Overall, our results are quite similar to the corresponding experimental data for DNA translocation through a nanopore [24], which we discuss in more detail in the following section.
At a next step, our data were shifted and scaled so that the distribution curve starts at zero-time and the total probability is equal to unity. The resulting distributions are on average not gaussians, but skewed towards long translocation times, consistent with experiment [24]. Therefore, the translocation time for each length is not assigned to the mean time, but to the most probable time, which is the position of the maximum in the histogram. In Fig.4 the distribution of all the events for polymer sizes N=50, 100 and 300 are shown. In this figure, the most probable translocation time for each length is denoted by an arrow. From this analysis, a nonlinear relation between the most probable translocation time τ 0 and the polymer length is obtained that follows closely the theoretically expected scaling τ X (N ) ∼ N α , with α ∼ 1.29 (see Fig.4).  3.4. Dynamics with and without a solvent. A closer inspection into the polymer dynamics reveals some interesting features. The molecule shows a blob-like conformation on either side of the membrane as it moves through the hole. It may either translocate very fast or move from one chamber to the other intermittently, with pauses. Both types of events are present with and without a fluid solvent. In addition, a careful analysis of all the translocated chains unravels the difference between slower and faster translocation within the same fast translocation regime. The nature of the variations in time is connected to the random fluctuations of the polymer throughout its motion, rather than the temperature or its length. These fluctuations are correlated to the entropic forces (gradient of the free energy with respect to a conformational order parameter, typically the fraction of translocated beads, see r(t) below) acting on both translocated and untranslocated parts of the polymer. In fact, when a solvent is present, the interplay between these forces and F drag , F pull determines the motion and the shape of the chain and thereby the translocation time. At some point part of the chain shapes up in an almost linear conformation increasing in this way the entropic force acting on it. This eventually leads to deceleration of the whole chain. Fig.6 shows an illustration of this argument, where a polymer chain, surrounded by a solvent is represented at a time where it starts to slow down. In this figure, a polymer with the same length but different initial configuration is also shown at the same time. It is very instructive to monitor the progress in time of the number of translocated monomers N (t). Note that r(t) ≡ N (t)/N serves as a reaction-coordinate, with the translocation time defined by the condition r(t X ) = 1. The translocated monomers for processes with and without hydrodynamics are shown in Fig.7. For the former, events related to the polymers of Fig.6 are shown (curves A 1 , A 3 ), as well as that related to the most probable time (A 2 ). The arrow in this figure indicates the timestep corresponding to the snapshots in Fig.6. The translocation for a given polymer proceeds along a curve virtually related to its initial configuration and its interactions with the fluid. It is clearly visible that there is no general trend. The non-hydrodynamic case is in principle different, especially in terms of the time range which is larger. This reveals the importance of hydrodynamic coherence.
Additional insight into the dynamics is obtained by altering the parameter set. This has not yet been extensively explored, but it was found that a choice of k B T = 10 −5 , F pull = 0.01, and ǫ = 0.002 leads to the frequent retraction of the polymer. In other words, after having translocated a large fraction of its length, the polymer occasionally reverses its motion and anti-translocates away from the hole, never to find its way back into it. Moreover, we find that a polymer that retracted in the presence of a solvent, manages to fully translocate if the solvent is absent. It is interesting to observe that, in principle, no such type of anti-translocating behavior has been observed for short polymers. This indicates that hydrodynamics significantly speed-up and alter the nature of translocation, especially for long polymers at low temperatures. This highly irregular dynamics escapes any scaling or statistical analysis, as well as dynamic Monte Carlo simulations [25], and can only be revealed by self-consistent many-body hydro-dynamic simulations.
4. DNA translocation through a nanopore. The translocation of biopolymers, such as DNA and RNA plays a major role in many important biological processes, such as viral infection by phages, inter-bacterial DNA transduction, and gene therapy [26]. The importance of this process has spawned a number of in vitro experiments, aimed at exploring the translocation process through micro-fabricated chan- nels under the effects of an external electric field, or through protein channels across cellular membranes [27]. In particular, recent experimental work has focused on the possibility of fast DNA-sequencing by reading the base sequence as the polymer passes through a nanopore. Some universal features of DNA translocation can be analyzed by means of suitably simplified statistical schemes [28] and non-hydrodynamic coarsegrained or microscopic models [29,30]. However, a quantitative description of this complex phenomenon calls for state-of-the art modeling of the type described above. Accordingly, we explore here to what extent the results discussed above for the generic situation of polymer translocation apply to the DNA case.
First, we note that, as already mentioned in the previous section, our results are quite similar to the experimental data for DNA translocation through a nanopore [24]. Three different interpretations of the current model are physically plausible: (a) Following the framework used in recent studies of DNA packing in bacteriophages [31], one monomer in our simulation can be thought of as representing a DNA segment of about 8 base-pairs, that is, each bead has a diameter of 2.5 nm, the hydrated diameter of B-DNA in physiological conditions. (b) It is also physically plausible to assume that a bead represents a portion of DNA equivalent to its persistence length of about 50 nm, which translates into mapping one bead to ∼150 base-pairs. (c) Alternatively, as is typically done in simulations of the λ-phage DNA in solution [32], one bead can be taken to correspond to ∼ 10 3 base-pairs. In all three cases h = ∆x is equal to the bead size, while the pore, having a width of 2∆x, will be different from the pores used experimentally, either smaller or larger. In addition, the coarse graining model that handles the DNA molecules indicates that the MD timescale is stretched over the physical process. A direct comparison between our probability distributions for polymer translocation and the experimental results sets a different MD timestep for the cases (a), (b), and (c) which is of the order of 3 nsec, 100 nsec and 5 µsec, respectively, leading to a LB timestep ∆t = 5dt of 15 nsec, 500 nsec, and 25 µsec. It is difficult at this stage to assign a unique interpretation of our model in relation to a physical system. A thorough exploration of the parameter space is required before such an assignment can be made convincingly. This is beyond the scope of the present work but will be reported in future publications.
A second encouraging comparison is that the scaling we found for the translocation time with polymer length (with exponent α = 1.29) is quite close to the experimental measurement for DNA translocation [24] (α = 1.27 ± 0.03). Beyond the apparent consistency between experiment and theory, additional insight can be gained by analyzing the polymer dynamics during translocation.
A hydrodynamic picture of DNA translocation has been presented in ref. [24]. In this work, the authors assume that the electric field drive is in balance with the Stokes drag exerted by the solvent on the blob configuration of the polymer, that is where ρ is the density, ν the kinematic viscocity, τ the translocation time, and R U the translocated part of the radius of gyration. In order for this balance to apply at all times, it is clear that R 2 must be constant in time, hence it cannot be identified neither with the translocated nor with the untraslocated gyration radius of DNA. To this end, in Fig.8 we represent the time-evolution of the radii of gyration for the two sections of DNA, R U (t) and R T (t), the untranslocated (U ) [with x > hN x /2] and translocated (T ) [with x < hN x /2] parts, respectively. In order to identify a time-invariant radius, we define the R I (t) = QN ζ I (t), where I = U, T stands for the untranslocated and translocated parts and N I (t) is the corresponding number of monomers. The exponent ζ ∼ 0.6 is the same as previously noted and Q is a constant (for long enough polymers). If translocation could be described by the dynamics of a single-blob object, characterized by an effective radius of gyration, defined as then this quantity should be constant in time. Since the N = N T (t) + N U (t) holds for all t, the above relations lead to We focus first on the case without hydrodynamics, Fig.8(a). For very small chains R I (t) does not scale as N 0.6 I and the definition for R E is not valid at the first and last ∼ 15 − 20% parts of the event, during which the untranslocated and the translocated parts, respectively, are small. Outside these limits, R E (t) as obtained from the definition (4.1), with the values of R U (t), R T (t) directly taken from the simulations, is indeed approximately constant. In addition, the values R T (t = t X ) and R U (t = 0) do not coincide, since the former is lower than the latter. This is also the case when a solvent is present, Fig.8(b). Thus, regardless of the dynamic pathway and the different conformations the chain may possess during translocation, once the event is completed the polymer is more compact than at t = 0. Comparison of the cases with and without hydrodynamics reveals that in the latter case the polymer becomes up to ∼ 7% more confined than when a solvent is added. The untranslocated part of the radius of gyration at the end of the process shows an abrupt drop. As a consequence the polymer t = t X does not fully recover its initial volume. It it plausible, that by allowing the polymer to further advance in time, R T (t = t X ) will become similar to R U (t = t 0 ), but this remains to be examined. Nevertheless, in this work we have been interested mainly on the chain dynamics related to the first passage times, which correspond to the exact period of time needed until all the beads have translocated.

Conclusions.
We have presented a multiscale methodology based on the concurrent coupling of constrained molecular dynamics for the solute biopolymers with a lattice Boltzmann treatment of solvent dynamics. Owing to the dual field-particle nature of the Lattice Boltzmann technique, this coupling proceeds seamlessy in time and only requires standard interpolation/extrapolation for information-transfer in physical space. This multiscale methodology has been applied to the case of polymer translocation through a nanopore, with special emphasis on the role of hydrodynamic coherence on the dynamic and statistical properties of the translocation process. It is found that hydrodynamic interactions play a major role in accelerating the translocation process, especially for long molecules at low temperature.
An attempt to connect these results to the process of DNA translocation through a nanopore revealed certain similarities with experiment, especially in the scaling law of the translocation time with polymer length. The presence of hydrodynamic interactions lead to a decrease in the translocation times, compared to the cases without a fluid solvent. Inspection of the variation of the translocated beads and the radii of gyration with time reveals interesting aspects of the DNA dynamics during translocation.
Future directions for the simulations include the detailed study of the effects of temperature, finite-length and geometrical details of the nanopore geometry, as well as electrostatic interactions of the DNA molecule with the surrounding fluid. To this end, resort to parallel computing is mandatory, and we expect the favourable properties of LB towards parallel implementations to greatly facilitate the task. Work along these lines is currently in progress.