## Abstract

The effect of a slit-like confinement on the relaxation dynamics of DNA is studied via a mesoscale model in which a bead and spring model for the polymer is coupled to a particle-based Navier–Stokes solver (multi-particle collision dynamics). The confinement is found to affect the equilibrium stretch of the chain when the bulk gyration radius is comparable to or smaller than the channel height and our data are in agreement with the (*R*_{g,bulk}/*h*)^{1/4} scaling of the polymer extension in the wall tangential direction. Relaxation simulation at different confinements indicates that, while the overall behaviour of the relaxation dynamics is similar for low and strong confinements, a small, but significant, slowing of the far-equilibrium relaxation is found as the confinement increases.

## 1. Introduction

A great number of microfluidic and nanofluidic devices have been proposed over the last few years for separation and manipulation of large biomolecules. In this context, the understanding of the dynamics of a single DNA molecule in confined geometry deserves greater interest because of the possible implications for design and improvement of devices that implement biological processes.

After the recent development of appropriate techniques to visualize single macromolecules, a number of experiments analysing the effect of confinement on a single DNA molecule have appeared; namely, the equilibrium properties of confined DNA have been studied by several groups (e.g. [1–4]). Attention has been mainly devoted to the modification of certain observables such as the equilibrium stretch, diffusion constant and relaxation time, and their resulting scalings depending on the degree of confinement interpreted in the framework of blob theory [5]. In the case of slit-like geometry, if the channel height *h* is comparable to, or smaller than, some relevant equilibrium scale then confinement breaks isotropy and, as reported by Balducci *et al.* [3], other alterations of the dynamical behaviour set in. In particular, this recent study shows that the dynamics of a stretched DNA molecule in a slit microchannel is characterized by two distinct relaxation regimes, the slower one emerging only when the channel height is smaller in dimension than the tension blob.

Numerical investigation of confined DNA molecules has been recently carried out using a different approach. On the length scale of the confinement experiment, the liquid can be treated as a continuum [6,7]; however, to analyse the single polymer dynamics, a discrete, though coarse-grained, model for the macromolecule is required and both hydrodynamics and thermal coupling need to be accounted for.

Jendrejack *et al.* [8] proposed a bead and spring Brownian model in which hydrodynamics is accounted for via the Oseen tensor. However, this approach is difficult to apply in confined geometries because the hydrodynamics interaction tensor cannot be calculated analytically and a numerical procedure is required (see also [9,10]). To overcome this difficulty, many groups have used an appropriate beads–springs model with a coarse-grained solvent whose momentum transport is explicitly computed (e.g. [11–13]). In the latter, the bead-and-spring model proposed by Jendrejack *et al.* [8] is coupled with the multi-particle collision dynamics (MPCD) approach for the solvent developed by Malevanets & Kapral [14]; see also Gompper *et al.* [15] for a thorough account of the applications. This method (sometimes called stochastic rotation dynamics) has been used for the last few years to analyse single polymer molecules in a solvent. In particular, equilibrium properties in non-confined geometries have been verified with a careful investigation of the model parameters by Mussawisade *et al.* [16] and Ripoll *et al.* [17].

In this paper, we analyse the confined equilibrium and non-equilibrium properties of a single DNA molecule. The bulk quantities (used as references for the confined cases) have been obtained via unconfined simulations not shown here. The polymer extension is found to significantly feel the slit-like confinement only when the bulk gyration radius is equal to or smaller than the channel height. Relaxation simulation was performed starting with an elongated configuration. It was found that the overall behaviour of the relaxation curves is similar for low and high confinement; however, a significant slowing of the far-equilibrium relaxation dynamics is apparent as the confinement increases.

## 2. The model

To construct models for these systems, it is crucial to coarse-grain the details and to keep the relevant length and time scales in order to observe the phenomena in which we are interested. The validity of the present approach has already been tested by other authors [13,16,17]; however, we will briefly recall the details of the numerical technique in the present section in order to highlight the differences from previous studies.

### (a) Solvent dynamics

For modelling the solvent, the MPCD approach developed by Malevanets & Kapral [14] is used. In MPCD, the solvent is modelled as a system of *N* particles of mass *m*, position **x**_{i} and velocity **v**_{i}, *i*=1,…,*N*. Time is discretized in intervals of length Δ*t*. The evolution of the system through a single time interval Δ*t* is the result of two steps, namely streaming and collision. In the streaming step, the particle dynamics is governed by Newton’s Second Law. The solvent particles do not interact with each other; hence, in the absence of external forcing, the evolution is simply
2.1
In the collision step, the system is divided into cubic cells of side *l*. For each cell c, the mean velocity **V**_{c} is calculated as the average of instantaneous velocities of the particles inside the cell. Each particle velocity is then changed according to where the superscripts a and b stand for ‘after’ and ‘before’ collision, while *ω* is a random rotation matrix. Here, we use a constant rotation angle, *α*=5/6*π*, in a uniformly distributed random direction. Several theoretical and numerical studies have analysed the properties of MPCD solvents. In particular, MPCD is able to reproduce on average the hydrodynamic motion [14,18], and the transport coefficients, such as viscosity, can be estimated as a function of statistical properties of the random matrix and fluid properties (density, temperature and particle mass) as reported by Kikuchi *et al.* [19] and Ihle *et al.* [20]. An additional step is implemented to ensure Galilean invariance [21]; namely, the grid defining the collision cells is shifted before each collision by a random vector whose components are extracted from a uniform distribution in the [−*l*/2,*l*/2] interval. In the present study, we use the length of the MPCD cell *l* as the length unit, the mass *m* of a solvent particle as the mass unit and *k*_{b}*θ* as the energy unit. The collision time step is Δ*t*=0.1, with a resulting viscosity (obtained from the expression in [19]) of *μ*=10.0042.

### (b) Polymer dynamics

The polymeric chain is divided into segments, each representing a subset of monomers and modelled by a bead and a spring. With such a choice, the description of physical short-range interactions is neglected to achieve a feasible model for long-range motions. In detail, our polymer is modelled, following Jendrejack *et al.* [8], as a sequence of *N*_{b} beads of mass *M* whose positions are indicated as **R**_{i}. Adjacent beads interact by means of the spring potential
2.2
where *b*_{k} is the Kuhn length and **R**_{ij}=**R**_{i}−**R**_{j}. *N*_{k} is the total number of Kuhn segments, *N*_{k,s} is the number of segments in each spring and *N*_{s}=*N*_{b}−1 is the number of springs; *R*_{0}=*N*_{k,s}*b*_{k} and the total contour length of the polymer *L*_{c}=*N*_{k}*b*_{k}=*R*_{0}*N*_{s}. To appropriately model a polymer in a dilute solution both the excluded volume and the hydrodynamic interaction (HI) must be taken into account. The latter is implicitly embedded in this approach by the presence of the solvent with the choice of parameters used, while the excluded volume interaction between all the bead pairs is
2.3
with
2.4
and *v* the excluded volume parameter. In the aforementioned units, the mass of a single bead is *M*=10, while the other parameters are *b*_{k}=0.14, *N*_{k,s}=17.4 and *v*=0.012. The beads do participate in the collision phase; in particular, the cell mean velocity **V**_{c} is calculated as the velocity of the centre of mass of the cell. These parameters are taken from Watari *et al.* [13], who found a reasonable agreement with good solvent Zimm scaling for the gyration radius, diffusion coefficient and relaxation time in the bulk, for polymer chains of *N*_{b}=5÷80. Equations of motion for the beads are numerically integrated using a standard velocity Verlet scheme with time step, *δt*=0.1Δ*t*.

### (c) Confined simulations set-up

We performed simulations in slit-like geometry with the two walls parallel to the *Oxy* plane; the *z*-coordinate of the wall is indicated as *z*_{w}. In the streaming phase, the standard bounce-back rule is applied to the solvent particles. The polymer confinement is obtained with the repulsive Lennard–Jones (LJ) potential (*σ*_{w}=1 and *ϵ*_{w}=1)
2.5
for *z*_{i}−*z*′_{w}<2^{1/6}*σ*_{w} and zero for (*z*_{i}−*z*′_{w})>2^{1/6}*σ*_{w} acting on each bead. In expression (2.5), *z*_{i} is the coordinate of the bead and *z*′_{w} is that of a plane placed at a distance *σ*_{w} inside the wall position (i.e. the *z*=*z*_{w} for the plane at which bounce-back is performed for solvent particles). Actually, with this choice, a bead at the wall (*z*_{i}=*z*_{w}) experiences a potential of *k*_{b}*θ*; hence, owing to the stiffness of LJ potential, the beads hardly cross the wall. On the other hand, the LJ potential affects only a small slab (2^{1/6}−1)*σ*_{w}≃0.12 thick, close to the wall, resulting in the fact that the polymer is able to explore the whole channel (figure 1). Concerning the collision phase, we follow the implementation of Lamura *et al.* [18], where virtual thermalized particles are added to the collision cells intersected by the wall.

## 3. Results

To quantitatively study the effect of confinement on relevant statistical observables, a set of simulations were carried out in slit-like geometries for a polymer with *N*_{b}=80. We report the results obtained for four different channel heights, namely *h*=18,9,7,5. In figure 2*a*, the average values of the equilibrium extensions in the wall normal direction (*Z*_{ex,c}) and in the wall parallel direction (*W*_{eq,c}=(*X*_{eq,c}+*Y* _{eq,c})/2), normalized with their respective free space values, are reported as functions of the ratio between the channel height and the unconfined equilibrium extension, *h*/*X*_{ex,bulk}, whose value is roughly half of the dimensionless channel, *h*/*R*_{g,bulk}. It is apparent that, for *h*<*X*_{ex,bulk}, the extension in the *z* direction is smaller than that of the free space, and the one parallel to the wall is slightly larger. For *h*>2*X*_{ex,bulk}, the isotropy is fully recovered. This result is in agreement with the findings reported in the literature; see, for example, the study of Jendrejack *et al.* [9] and Chen *et al.* [2], in which the same transition as in figure 2*b* between a so-called weak confinement and a strong one has been observed when the inverse of the dimensionless channel (i.e. *R*_{g,bulk}/*h*) is larger than 0.4 for both microchannels and slit confinements.

According to Brochard & De Gennes [5], the behaviour of a linear polymer of persistence length *b*_{k} and contour length *N*_{k}*b*_{k} confined in a slit-like channel with height *h* between the equilibrium length and *b*_{k} can be described as follows. The molecule is seen as a chain of blobs of dimension *h*, whose configuration follows a two-dimensional self-avoiding walk and no HI is assumed. The segments of the chain inside the blob follow a three-dimensional self-avoiding walk and they interact hydrodynamically. Under these assumptions, we obtain a scaling for the equilibrium extension with respect to height *h* as
3.1
where the values of the exponents have been obtained using the ideal values 3/4 and 3/5 for the Flory–Edwards exponents *ν*_{2D} and *ν*_{3D}, respectively [22]. As shown in figure 2*b*, in our set-up, the cross-over between weak and strong confinement occurs for two of the four simulations, and for those the comparison with the expected scaling regime (3.1) is good within the confidence interval. However, with the current choice of the simulation parameters, we do not reach a very high confinement regime, as in Chen *et al.* [2]. Further simulations with longer polymers have been planned. It is worth pointing out that, also from the derivation, the stretching scaling is only owing to excluded volume constraints, between the beads and with the wall.

The relaxation dynamics is studied via non-equilibrium runs in which the DNA molecule, initially in elongated conformation along the *x* direction, spontaneously coils up until it reaches an equilibrium configuration. Independent randomly generated configurations, in which the polymer elongation in the *x* direction was about 60 per cent of the polymer contour length *L*_{c}, were used as initial conditions for the runs. The ensemble averages estimated on at least 20 independent simulations for each channel height are calculated for the time evolution of polymer extension.

In order to highlight the characteristics of the relaxation dynamics better at different channel heights, following Balducci *et al.* [3], we analyse the average scaled-squared extension along the *x* direction defined as . Figure 3 reports the time evolution of the scaled extension for *N*_{b}=80 and *h*=5,7,9,18 where the time is set to zero when the polymer extension is 30 per cent of the contour length (〈*X*_{ex}(*t*)/*L*_{c}〉=0.3). The overall behaviour of the curves is very similar; however, a slight difference in the far-equilibrium relaxation is observed (inset of figure 3). In particular, the narrower the channel, the slower the relaxation; hence, the confinement appears to slow the relaxation process. This observation confirms that HI is implicitly taken into account in our model and hence the presence of the wall modifies the dynamics even in the early stages of relaxation. For a complete discussion see Jendrejack *et al.* [9], in which the comparison of a free-draining model and that with explicit hydrodynamic interaction allows us to attribute the modification of the relaxation dynamics to the confinement only when HI is considered. In figure 3*a*, the more confined case is *h*=5. To analyse the modified relaxation dynamics far from equilibrium and its dependence on the confinement *h*, we have rescaled the abscissa of the other three relaxation curves by a factor *τ* in order to collapse on the *N*_{b}=80, *h*=18 curve. Such a value is a measure of alteration on the relaxation dynamics. The scaling of such a parameter with the degree of confinement is reported in figure 3*b*, where the scaling exponent is extracted and is equal to 0.42. It is worth pointing out that this value is larger than that obtained by Jendrejack *et al.* [9], which was calculated under near-equilibrium conditions, and smaller than that from Balducci *et al.* [3], for the far from equilibrium conditions. These two other scalings are also reported in figure 3*b* as a reference. Namely, Balducci *et al.* [3] in a slit-like experimental set-up, which is qualitatively reproduced by the current simulations, were able to identify two different scaling regimes at larger extensions and near equilibrium, respectively. Moreover, the authors suggest that the direct observation of the second regime is possible only if the confinement is large enough that the equilibrium stretch 〈*X*^{2}_{eq,c}〉 is smaller than a predicted cross-over length. In the present simulations with *N*_{b}=80, we do not meet this criterion.

## 4. Concluding remarks

In this paper, the effect of a slit-like confinement on both equilibrium and non-equilibrium dynamics of a DNA molecule is studied. The results obtained for the far from equilibrium regime of the relaxation dynamics are similar to the scaling observed in Balducci *et al.* [3]. The predictive capability of the present approach suggests that further investigations are advisable to address the dynamics of the relaxation in near-equilibrium conditions and to test the results against blob theory predictions.

## Acknowledgements

Computing resources were made available by CASPUR under HPC grant 2009.

## Footnotes

One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.

- This journal is © 2011 The Royal Society