BROWNIAN DYNAMICS SIMULATION OF DNA IN COMPLEX GEOMETRIES

by Yu Zhang

A dissertation submitted in partial fulfillment of the requirements for the degree of

Doctor of Philosophy (Chemical Engineering)

at the UNIVERSITY OF WISCONSIN – MADISON 2011

i

To my wife, Xi For your love, your support, and a wonderful child.

To my son, Alexander For your laughters, which brighten my life.

To my parents, Liangzhu Zhang and Zhenglan Yu For your unconditional love.

ii

Acknowledgments Thanks to my advisors, Michael Graham and Juan de Pablo for their support, patience, and guidance. Thanks to former colleague Hongbo Ma, Juan Hernandez-Ortiz, and Patrick Underhill for many useful discussions on fluid mechanics and polymer physics. Thanks to current colleague Amit Kumar and Pieter Janssen for their helps and many discussions. Thanks also to other former and current colleagues of the Graham group – Wei Li, Li Xi, Samartha Anekal, Matthias Rink, Mauricio Lopez, Aslin Izmitli, Pratik Pranay, Kushal Sinha, Rafael Henriquez Rivera, and Friedemann Hahn – for all your helps and many discussions. Research projects presented in this dissertation are financially supported by the University of Wisconsin-Madison Nanoscale Science and Engineering Center funded by National Science Foundation.

iii

Abstract This dissertation is concerned with the dynamics of a long DNA molecule in complex geometries, driven by either electrostatic field or flow field. This is accomplished primarily through the use of Brownian dynamics simulation, which captures the essential physics at mesoscopic length scale, and allows us to simulate events happening on long time scale, such as DNA pore translocation and cyclic dynamics of a tethered DNA molecule in shear flow. Our work is novel in that, accurate and efficient algorithms have been designed and developed for both electric field-driven and flow-driven systems in complex geometries. We have focused on three distinct problems, which we describe below. In the electric field-driven case, we propose a novel class of electric field-actuated soft mechanical control element for microfluidics. This type of element employs the idea that under confinement a single polymer molecule is essentially a nanoscale porous media. The chain could block the passageway of relatively large analytes such as cells. At the same time, polyelectrolyte molecules, such as DNA, could deform and squeeze through narrow pores when a sufficiently strong electric field is applied. Brownian dynamics (BD)/Finite Element Method (FEM) simulation efficiently explore the design space, and results demonstrate that the On/Off switching could be

iv achieved within a proper parameter space. To study the effects of a solid impenetrable wall on the dynamics of a nearby DNA molecule, we examine the cyclic dynamics of a single DNA molecule tethered to a hard wall in shear flow, using Brownian dynamics simulation. We focus on the dynamics of the free end (last bead) of the tethered chain and we examined the crosscorrelation function and power spectral density of the chain extensions in the flow and gradient directions as a function of chain length N and dimensionless shear rate W i. Extensive simulation results suggest a classical fluctuation-dissipation stochastic process and question the existence of periodicity of the cyclic dynamics, as previously claimed. We support our numerical findings with a simple analytical calculation for a harmonic dimer in shear flow. In the case of flow-driven DNA molecule in micro/nano-fluidics, one big challenge is that an efficient algorithm is required to calculate fluctuating hydrodynamic interactions (HI) in complex geometries. We have developed an accelerated immersed boundary method that allows fast calculation of Brownian motion of polymer chains and other particles in complex geometries with HI. With this new method, the first detailed analysis of a recent set of interesting nanofluidic experiments involving DNA dynamics in a complex flow geometry is performed. This analysis explains the observed dynamics over a wide range of parameter values (flow rate, molecular wieght) and illustrates the important quantitative effect of the hydrodynamic interactions on the behavior of the system.

v

Contents Acknowledgments Abstract List of Figures

ii iii viii

List of Tables

x

1 Introduction

1

2 Problem Statement

4

2.1

Bead-spring DNA Model . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2

Kinetic theory of single polymer molecule . . . . . . . . . . . . . . . .

6

2.3

The Langevin equation and Brownian dynamics simulation . . . . . .

9

3 Intramolecular interactions and external forces

12

3.1

Spring force law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

3.2

Excluded volume interactions . . . . . . . . . . . . . . . . . . . . . .

14

3.3

Hydrodynamic interactions . . . . . . . . . . . . . . . . . . . . . . . .

15

3.4

Electrophoretic force . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

vi 3.5

Polymer-wall steric interaction . . . . . . . . . . . . . . . . . . . . . .

4 Bistability and field-driven dynamics of confined tethered DNA

19 21

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

4.2

Polymer model and simulation approach . . . . . . . . . . . . . . . .

28

4.3

Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . .

31

4.3.1

Pore-crossing geometry . . . . . . . . . . . . . . . . . . . . . .

31

4.3.2

Pore-entry geometry . . . . . . . . . . . . . . . . . . . . . . .

37

Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

4.4

5 Tethered DNA dynamics in shear flow

43

5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

5.2

Discussion of Methods for Hydrodynamics of Polymer Solutions . . .

49

5.2.1

Implicit Solvent: Brownian Dynamics . . . . . . . . . . . . . .

50

5.2.2

Explicit Solvent: Continuum Methods . . . . . . . . . . . . .

52

5.2.3

Explicit Solvent: Particle Methods . . . . . . . . . . . . . . .

55

5.2.4

Coupled Methods . . . . . . . . . . . . . . . . . . . . . . . . .

56

Simulation Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

5.3.1

Brownian Dynamics . . . . . . . . . . . . . . . . . . . . . . .

57

5.3.2

Lattice-Boltzmann . . . . . . . . . . . . . . . . . . . . . . . .

60

5.3.3

Stochastic Event-Driven Molecular Dynamics . . . . . . . . .

62

5.4

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

5.5

Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

5.3

6 Flow-driven DNA dynamics in complex geometries 6.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80 81

vii 6.2

Methods for hydrodynamics of confined polymer solutions . . . . . .

86

6.3

Polymer model and simulation method . . . . . . . . . . . . . . . . .

91

6.3.1

Model and governing equations . . . . . . . . . . . . . . . . .

91

6.3.2

Governing equations . . . . . . . . . . . . . . . . . . . . . . .

92

6.3.3

Mobility tensor and Fixman’s midpoint algorithm . . . . . . .

94

6.3.4

Chebyshev approximation . . . . . . . . . . . . . . . . . . . .

96

6.3.5

Fast Stokes solver with complex boundary conditions . . . . .

97

6.4

6.5

DNA flowing across an array of nanopits . . . . . . . . . . . . . . . . 114 6.4.1

Dynamics at low Peclet number . . . . . . . . . . . . . . . . . 117

6.4.2

Dynamics at high Peclet number . . . . . . . . . . . . . . . . 123

Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

7 Conclusion and future work

126

A Lattice random walk model of a tethered polymer

129

B BD/FEM algorithm for simulating DNA electrophoresis

137

C Dimer in shear flow

141

D Fixman’s midpoint algorithm

146

viii

List of Figures 2.1

Schematic of a bead-spring chain in Brownian dynamics simulation. .

6

4.1

Schematic representations of soft nanomechanical bistable elements. .

26

4.2

Soft nanomechanical elements. . . . . . . . . . . . . . . . . . . . . . .

28

4.3

Top view of the simulation domains. . . . . . . . . . . . . . . . . . .

29

4.4

Time trajectories for a 3D pore-crossing system . . . . . . . . . . . .

32

4.5

Snapshots of a pore-crossing event. . . . . . . . . . . . . . . . . . . .

33

4.6

Simulation results for the quasi-2D and 3D pore crossing geometries.

35

4.7

Time trajectory and probability density for the pore entry geometry .

38

4.8

Simulation results for the quasi-2D pore-entry geometry . . . . . . . .

39

4.9

Phase transition of the “competitive bistability” system. . . . . . . .

40

5.1

Cyclic dynamics of a tethered DNA molecule . . . . . . . . . . . . . .

46

5.2

Probability distribution of the end bead of the tethered DNA molecule

68

5.3

Relaxation time of the tethered DNA . . . . . . . . . . . . . . . . . .

70

5.4

Spectral analysis of the end-bead position time series . . . . . . . . .

73

5.5

Cross correlation functions for different beads . . . . . . . . . . . . .

74

5.6

Comparison of the cross-correlation function among different methods

76

6.1

DNA flowing through nanopits . . . . . . . . . . . . . . . . . . . . . .

85

ix 6.2

Schematic of the immersed boundary method . . . . . . . . . . . . .

99

6.3

Properties of the mobility tensor for the boundary particles . . . . . . 103

6.4

Screening function for the GGEM algorithm . . . . . . . . . . . . . . 105

6.5

Comparison between delta and regularized delta point forces . . . . . 106

6.6

Comparison between Hasimoto’s solution and GGEM . . . . . . . . . 110

6.7

Error analysis of GGEM . . . . . . . . . . . . . . . . . . . . . . . . . 111

6.8

Laminar channel flow . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

6.9

Uniform flow around a sphere . . . . . . . . . . . . . . . . . . . . . . 115

6.10 Schematic of the immersed boundary representation of the device . . 116 6.11 Streamlines and contour plot of the streamwise velocity in a pit . . . 117 6.12 Snapshots of a hopping event . . . . . . . . . . . . . . . . . . . . . . 118 6.13 DNA dynamics at low Peclet number . . . . . . . . . . . . . . . . . . 120 6.14 Mean resident time τ v.s. Peclet number P e. . . . . . . . . . . . . . . 121 6.15 Mean residence time τ v.s. chain length N and Peclet number P e . . 122 6.16 DNA dynamics at high Peclet number . . . . . . . . . . . . . . . . . 124 A.1 A tethered polymer as lattice random walk. . . . . . . . . . . . . . . 133 A.2 The repulsive wall-polymer potential . . . . . . . . . . . . . . . . . . 135 B.1 Schematic of the nearest neighbor search algorithm. . . . . . . . . . . 140 C.1 Spectral analysis of a dimer in shear flow . . . . . . . . . . . . . . . . 145

x

List of Tables 2.1

Properties of λ-DNA.

. . . . . . . . . . . . . . . . . . . . . . . . . .

3.1

Spring models and the corresponding force laws.

. . . . . . . . . . .

5 14

1

Chapter 1 Introduction Transport of polymer solutions in constricted spaces [1] is a long-standing research topic with many applications including polymer enhanced-oil-recovery, size exclusion chromatography [153], gel-electrophoresis [151], and recently single DNA molecule analysis using micro- and nano-fluidic devices [140, 125]. In many fluidic devices designed for separation, manipulation, and sequencing of DNA, the critical dimension of the constriction approaches the size of the polymer molecule or smaller, and understanding the effects of spatial restrictions is essential to the rational design and applications of those devices. In particular, capturing the interplays between polymer dynamics, spatial confinement, and electrostatic or flow field is essential to a proper description of these devices. It is therefore of great interest to develop validated fast modeling and simulation tools to understand the fundamental physics of polymer solutions in complex geometries. There are two major types of driving forces to transport DNA molecules through microfabricated devices: electrokinetic forces (mainly electrophoretic force) and fluid

2 flow drag force. Many devices with well-defined microstructures have been proposed for DNA manipulation, especially for separation and stretching purposes using electrophoresis [45]. By contrast, much less attention has been paid on pressure (flow) driven DNA dynamics in microfabricated devices [141, 36, 143]. There are significant differences between the electrophoresis case and the pressure driven case. Consider for example DNA through a slit geometry. In electrophoresis, in a uniform field, the velocity of each DNA segment is the same everywhere in the channel; by contrast, the unperturbed velocity profile is parabolic across the channel in the pressure driven case at low Reynolds number, and this velocity gradient can give rise to interesting transport phenomena, such as Taylor dispersion [141]. Furthermore, hydrodynamic interactions (HI) between objects such as polymer segments in an unconfined domain are long-ranged, leading to strong many-body effects, for example, Zimm scaling of the self-diffusion coefficient of polymer. In confined geometries, the long-ranged nature of hydrodynamic interactions changes substantially, leading to significant changes in polymer dynamics. [1]. In a slit geometry, for example, the velocity field due to a point force perpendicular to the walls decays exponentially. For a force parallel to the walls, the velocity field has a parabolic form in the wall-normal direction and decays as 1/r 2. This is still a slow decay, but the symmetry of the flow leads to cancellations upon averaging that result in screening [9, 144, 14]. Both experiments and simulations of the diffusion of long flexible DNA molecules in slits [31, 78, 75] are consistent with screening of hydrodynamic interactions on the scale of slit height. Nevertheless, on smaller scales, HI are not screened and lead to changes in segment mobilities and boundary effects such as cross-stream migration, which leads to depletion layers much larger than the equilibrium chain size [77, 90, 30, 68, 29, 81]. A

3 significant portion of this thesis is concerned with developing general methodologies for Brownian dynamic simulations of DNA in complex geometries, driven by either electrostatic field (Chapter 4 ) or flow field (Chapter 5 and Chapter 6). This dissertation is organized as follows: In Chapter 2, we present the governing equations for the dynamics of a single DNA molecule. Intramolecular interactions and external potentials are introduced in Chapter 3. In Chapter 4 we present the algorithm for DNA electrophoresis in complex geometries as applied to study the properties of a class of soft nanomechanical control elements for microfluidics. In Chapter 5, we discuss the Brownian dynamics simulation results of the cyclic motion of a DNA molecule which is tethered onto a solid wall and experiences a shear flow field. In Chapter 6, we develop and validate a method for the efficient calculation of fluctuating hydrodynamics in complex geometries. In that chapter we also present simulation results for the dynamics of a flow-driven DNA through a nanofluidic slit with an embedded array of nanopits. Conclusions and future work are given in Chapter 7. The different sections in Chapters 4 through 6 correspond to different publications and manuscript. Those sections are therefore fairly self-contained, and some repetition should be expected.

4

Chapter 2 Problem Statement In this dissertation, we are concerned with the dynamics of a single DNA molecule which is confined in a fluid domain with complex boundary conditions. As we are interested in long-time ( > 1 second) dynamics of the DNA molecule, a bead-spring DNA model and Brownian Dynamics simulation are used in this work. A general discussion of the wide range of techniques for modeling the hydrodynamics of polymer chains in solution is given in Chapter 5, and a brief overview of various methodologies for modeling hydrodynamics of polymer solutions in complex geometries is given in Chapter 6. In the remainder of this chapter, we describe the DNA model and the governing equations.

2.1

Bead-spring DNA Model

Double-stranded DNA is a semiflexible polymer whose physical and chemical properties have been studied extensively. The physical properties of λ-DNA, the bead-spring model of which is used in this work, are summarized in Table 2.1. The choice of DNA

5 Table 2.1: Properties of λ-DNA. Property

Symbol

Value

References

Number of base pares

Nbp

48502

[137]

Contour length

Lc

22 µm

[137]

Persistence length

lp

53 nm

[138]

Kuhn length

bk

2lp

Radius of gyration

Rg

730 nm

[137]

Diffusion coefficient

D

0.47 µm2 /s

[137]

Electrophoretic mobility

µ0

4.2 × 10−4 cm2 /Vs

[142, 102]1

model depends upon the molecule properties (chemical, mechanical, etc.) one wants to model and the level of molecular detail one needs to retain [2]. For long time Brownian dynamics simulation, one of the most commonly used coarse-grained models for linear DNA molecules is the bead-spring model. The molecule is modeled as a series of Nb beads connected by Ns = Nb − 1 springs, as shown in Figure 2.1. The total number of degree of freedom for a freely moving bead-spring chain in therefore 3Nb . Bead-spring model is the most coarse grained version of a DNA model. The beads act as sources of fluid drag friction and the springs, which obey certain type of spring force law, represent a collection of persistence lengths. Carefully calibrated bead-spring chain model greatly enlarges the accessible time and length scales in numerical simulation. In this work, we use the model developed by Jendrejack et al.[76], which will be discussed in Chapter 3. 1

Measured in TBE buffer (0.01M TBE, pH = 8.3, at 23◦ Celsius).

6 i

3 ri 1

O

2

i+1 Nb

Nb-1

Figure 2.1: Schematic of a bead-spring chain in Brownian dynamics simulation. There are two equivalent ways to describe the dynamics of a DNA or, more generally, a polymer molecule. One is based on the Fokker-Planck equation (“Diffusion equation”) for the time evolution of the Nb beads phase space distribution function Ψ(R, t). The other approach uses Langevin equation to describe the motion of Nb beads. By solving the diffusion equation, one can obtain Ψ(R, t) directly. On the other hand, with the set of coupled Langevin equations, one can directly calculate trajectories of the beads and then obtain Ψ(R, t) by averaging. We will discuss first the diffusion equation, and then the Langevin equation and Brownian dynamics simulation.

2.2

Kinetic theory of single polymer molecule

Kinetic theory of polymer molecule describes the dynamics of a distribution function, i.e., evolution of the distribution function, under the action of forces. Using the bead-spring model of polymer, the positions of the beads represent the configuration of the polymer, and the configuration distribution function Ψ is defined so that ΨdR

7 is the probability to find the position vector of the first bead within dr1 around r1 , second bead at dr2 around r2 , and so on. The general equation for the configuration distribution function Ψ of a bead-spring polymer is called the Fokker-Planck equation or “Diffusion equation”. It is the derivation of this equation that we devote the remainder of this section. When the bead inertial relaxation times are short compared to the timescale of interest, it is often possible to ignore inertia, that is, we can write a force balance about bead i. The relevant forces on the bead i include: hydrodynamic force fih , Brownian force fib , and other non-hydrodynamic and non-Brownian forces fin which is the sum of spring force, excluded volume force, and steric force between the bead and wall. Then the force balance is written as fih + fib + fin = 0,

(2.1)

p fih = −ζ · [ui − (u∞ i + ui )],

(2.2)

in which

fib = −kB T fin = −

∂ ln Ψ, ∂ri

(2.3)

∂ φ. ∂ri

(2.4)

Equation 2.2 is the Stokes drag law: the hydrodynamic force on bead i is assumed to be linear in the slip velocity between the bead velocity ui and the solvent velocity p ∞ at the bead center (u∞ i + ui ) with the coefficient ζ. Here ui is the unperturbed

velocity when the polymer is absent, and upi is the perturbed velocity due resulting from the motion of other beads in the system. For now, we simply state that the

8 perturbed velocity at ri is linearly depend on the non-hydrodynamic, non-Brownian P b n forces acting on all beads in the system fin as upi = N j=1 Mij · fj . Here Mij is a tensor describing the hydrodynamic interaction between bead i and j.

Equation 2.3 gives the Brownian force. This equation has been derived by Bird et al. (Ref [21]) for the case of a structureless mass point in which the chain has been equilibrated in momentum space. Equation 2.4 represents the force from all the non-hydrodynamic and non-Brownian forces which stem from intra- and intermolecular potentials, and external potentials such as electric potential and bead-wall steric repulsive potential. We will return to a discussion of these potentials and forces in Chapter 3. Now substitute the expressions for various forces into the force balance equation and we obtain p − ζi · [ui − (u∞ i + ui )] − kB T

∂ ln Ψ + fin = 0. ∂ri

(2.5)

Rearranging above expression, we obtain 1 ∂ (−kB T ln Ψ + fin ) ζi ∂ri X 1 ∂ ln Ψ + fjn ). = u∞ ( δij I + Mij ) · (−kB T i + ζ ∂r j j j

p ui = (u∞ i + ui ) +

(2.6) (2.7)

Here, I is an identity matrix. From the first to the second line in above expression, we use the expression for upi . According to the kinetic theory, the configuration

9 distribution function satisfies the continuity equation X ∂ ∂Ψ =− ( · ui Ψ). ∂t ∂ri i

(2.8)

Defining the diffusion tensor D as

Dij =

kB T δij I + kB T Mij , ζj

(2.9)

and combining 2.7 and 2.8, we obtain the diffusion equation. X ∂ 1 X ∂ ∂Ψ · (u∞ [Dij · (−kB T ln Ψ + fjn )]Ψ)) =− ( i + ∂t ∂ri kB T j ∂rj i

2.3

(2.10)

The Langevin equation and Brownian dynamics simulation

In this section, we discuss the alternative descriptions of the dynamics of a polymer molecule in terms of random variable dW and its defining stochastic differential equation. This description is used in the Brownian dynamics simulation of a polymer molecule. The diffusion equation (Eq. 2.10) can be recast in the form of a Langevin equation (a stochastic differential equation)[105] dR = (U∞ + M · F + kB T B · BT = kB T M.

√ ∂ · M)dt + 2B · dW, ∂R

(2.11) (2.12)

10 Here R is the vector containing all bead positions {ri }, and the vector containing total non-Brownian, non-hydrodynamic forces acting on the beads is denoted F. The tensor B gives the magnitude of the Brownian displacement of the polymer beads, and is coupled to M by the Fluctuation-Dissipation theorem (Eq. 2.12). The vector dW is a random vector composed of independent and identically distributed random variables according to a real-valued Gaussian distribution with mean zero and variance dt. In Brownian dynamics simulation, applying a stochastic integration scheme to above equation generates the trajectories of a polymer molecule under forces. For example, using a simple forward Euler integration scheme we obtain R(t + ∆t) = R(t) + (U∞ + M · F + kB T

√ ∂ · M)∆t + 2B · ∆W. ∂R

(2.13)

In electrophoresis, DNA is a “free-draining” (FD) polymer. Therefore, we can ignore the hydrodynamic interactions (HI) between different chain segments, and the mobility tensor M is reduced to the product of

1 ζ

and an identity matrix. Hence, Eq. 2.13

is reduced to 1 R(t + ∆t) = R(t) + F∆t + ζ

s

2

kB T · ∆W. ζ

(2.14)

Above equation is relatively easy to solve compared with that in the flow-driven case. An efficient algorithm to find the electric field at the bead position in a field with steep gradient is essential. This issue is addressed in Appendix B. In the flow-driven case, to capture the full fluctuating hydrodynamic in complex geometries, we first need to turn Eq. 2.13 into a derivative-free form, as the expression for M is lacking when the molecule is confined in complex geometries. Then, the

11 stochastic differential equation can be solved based on a fast solver for Stokes flow in complex geometries. This will be discussed in Chapter 6.

12

Chapter 3 Intramolecular interactions and external forces The dynamics of the DNA molecule obeys Newton’s second law, and therefore, it is determined by the total force exerted on the chain. As discussed in last Chapter, the total force on a bead is composed of a viscous drag force f h , a Brownian force f b due to the collisions of the solvent molecules, and all other non-hydrodynamic nonBrownian forces f n , which includes spring forces, excluded volume forces, and any external forces such as electrophoretic force and DNA-wall steric interaction. In this chapter, we discuss those forces and potentials.

13

3.1

Spring force law

Each spring models a fraction of the full DNA molecule which is long enough such that its force-extension behavior obeys the Marko-Siggia spring law [26, 94], fs =

|re | −2 4|re | re kB T [(1 − ) −1+ ] , 2bk Nk bk Nk bk |re |

(3.1)

where f s is the average force needed to get an end-to-end vector of re for each DNA segment. The Kuhn length bk is twice the value of the persistence length, and when a DNA molecule is divided into Nk = Lc /bk segments, each Kuhn segments can be thought of as if they are freely jointed with each other (no bending, rotational, or torsional potentials). This spring law was derived based on a worm-like chain (WLC) model in polymer physics which is commonly used to describe the behavior of semi-flexible polymers. In this model, the molecule is treated as a flexible smooth rod. The rod’s local direction (tangential vector) decorrelates at distance s along the curve according to exp(−s/lp ), where the decay length, lp , is called the persistence length of the chain. The stiffer the chain, the larger the persistence length. For DNA, the persistence length is approximately 50 nm. The WLC model fits DNA force-extension data very well up to 10 pN, or < 95% stretching ratio [26]. As long as the number of Kuhn lengths contained in each spring is larger than 10bk , the Marko-Siggia spring law can be used [149]. In this work, each spring represents 20bk . There are other spring force laws proposed for simulation of polymer molecules, which are summarized in Table 3.1, together with the WLC model.

14 Table 3.1: Spring models and the corresponding force laws. Model

Force expression f s = Hre =

Hookean

fs =

FENE

fs =

Inverse Langevin model

3kB T r Nk b2k e

Hre 1−(|re |/Nk bk )2

kB T bk

L−1 ( N|rkeb|k ) |rree |

L(x) = cothx − WLC

3.2

fs =

kB T 2bk

[(1 −

|re | −2 ) Nk bk

1 x

−1+

4|re | re ] Nk bk |re |

Excluded volume interactions

In polymer physics, excluded volume interaction is a type of long-ranged1 interaction referring to the idea that in any real polymer molecule, two monomers cannot occupy the same space. This effect plays a far more important role in polymer solution than it does in solution of small molecules [35]. We often idealize the chain to allow overlap of monomers, and call it an ideal chain. By contrast, a regular chain with excluded volume interactions is called a real chain. Ideal chain does not exist in reality, but it is used extensively because it allows rigorous mathematical treatment of many questions. More importantly, the real chain behaves like an ideal chain in some situations, such as when the molecule is in concentrated solutions, melts, and in theta condition. In Appendix A, we derive an effective repulsive potential between a polymer molecule and a hard wall using the ideal chain model. In Brownian dynamics simulation for λ-DNA, we use an exponential form of the 1

Note that excluded volume interaction is long-ranged along the chain backbone. The potential is actually short-ranged in space (Eq. 3.2).

15 repulsive excluded volume potential between bead i and bead j, when we consider a real chain. It is obtained by considering the energy penalty due to the overlap of two Gaussian coil [41], −3|ri − rj |2 1 3 3/2 2 ) exp( ), Uijev = νkB T Nk,s ( 2 4πSs2 4Ss2

(3.2)

where ν is the excluded volume parameter, and Ss2 = Nk,s b2k /6 is the mean square radius of gyration of an ideal chain containing Nk,s Kuhn segments of length bk . The resulting expression describing the force acting on bead i due to the the presence of bead j within some cut-off distance is then 2 fijev = νkB T Nk,s π(

−3|ri − rj |2 3 5/2 ) exp( )(ri − rj ), 4πSs2 4Ss2

(3.3)

Although this potential is not self-consistent (any deformation of the coil caused by the overlap has been ignored), it does provide the correct scaling relationships [76].

3.3

Hydrodynamic interactions

Beads immersed in a fluid generate flows as they move due to various forces, and similarly they move in response to fluid motion through the Stokes drag. In this work, we treat bead i of the bead-spring chain as a sphere of hydrodynamic radius a. Then the relationship between the bead velocity ui and the drag force it exerts on the fluid is given by the Stokes law fi = ζ(ui − u(ri )), where ζ is the bead friction coefficient ζ = 6πηa, u(ri ) is the fluid velocity at the bead position, and η is the fluid viscosity. (Note that the finite size of the bead only arises in the friction coefficient.)

16 Through these hydrodynamic interactions (HI), beads interact with each other and with the walls of the confining geometry. In kinetic theory, the velocity perturbation due to the polymer molecule, up , is generally taken to be due to a chain of point forces acting on the fluid, and obtained by solving the incompressible Stokes flow problem, − ∇p(x) + η∇2 up (x) = − ∇ · up (x) = 0,

X

fn δ(x − xn ),

(3.4) (3.5)

up (r) = 0, r ∈ ∂Ωb ,

(3.6)

where p is the pressure, up is the fluid velocity, η is the fluid viscosity, fn is the force exerted on the fluid at point xn , δ(x) is the three-dimensional delta function, and ∂Ωb is the boundary of the fluid domain. In the actual implementation, regularized point forces are used, as further discussed flow. The velocity perturbation at the position of bead i due to the bead j is represented in a Greens function form as

upi =

Nb X j=1

Mij · fjn

(3.7)

where Mij is called hydrodynamic interaction tensor. To illustrate the components of the mobility tensor, we consider a two-bead chain (a dumbbell) in an unbounded domain, neglecting Brownian effects for the moment. Here the force balance shows that the fluid velocities u(ri ) experienced by bead i are given by 





u∞ 1





1 I ζ

∞

 



G (r1 − r2 )   f1     u1   , · +  =  1 ∞ f I G (r − r ) u∞ u2 2 2 1 2 ζ

(3.8)

17 where G∞ (r) =

1 (I 8πηr

+ rr/r 2) is the Oseen-Burgers tensor, the free-space Green’s

function for Stokes’ equations, r = |r|, and I is the identity matrix. The quantity G∞ (r1 − r2 ) · f2 gives the velocity at position r1 generated by the force f2 exerted by the particle at r2 on the fluid. We can write above expression succinctly as, U = U∞ + M · F. We note that

∂ ∂r

(3.9)

· M = 0 for G∞ . In simulation of coarse-grained polymer chain, the

singularity of the Oseen tensor must be regularized. The most common choice is the RPY tensor [120, 157]. In a confined geometry, a wall correction GW (r1 , r2 )has to be added to the mobility tensor M, and it becomes 

 M=

1 I ζ

+ GW (r1 , r2 )

∞



G (r1 − r2 ) + GW (r1 , r2 )  . 1 ∞ I + GW (r1 , r2 ) G (r2 − r1 ) + GW (r1 , r2 ) ζ

(3.10)

In general geometries, an analytical expression for GW (r1 , r2 ) is not available, but has to be calculated numerically. We note that it is not M that is needed, but rather, the product M · f which is the velocity generated by the total non-Brownian, nonhydrodynamic forces acting on the beads. As a result, using an accelerated Immersed Boundary Method, we can calculate M·f for polymer in complex geometries efficiently. We will return to a discussion on it in Chapter 6. In Green’s function-based representation of the hydrodynamic interactions, we assume that the beads are coupled instantaneously. In other words, the characteristic time for the hydrodynamic interactions to propagate through the system should be

18 much smaller than the time scale of other signals. These conditions are equivalent to that the Reynolds number is small Re ≪ 1 and the March number is small Ma ≪ 1, which are satisfied explicitly through the use of the Stokes equation.

3.4

Electrophoretic force

Electrophoretic force is commonly used to transport and manipulate DNA molecules. In this section, we discuss DNA electrophoresis and the effective charge assigned to each bead calculated using the electrophoretic mobility of a long DNA molecule. The velocity of a DNA molecule (> 100 bp) in free electrophoresis has little dependence on length [142]. When an electric field is applied across an electrolyte solution containing DNA molecules, the counter-ion cloud around the DNA experiences a force in the opposite direction to that of the DNA. This force generates ab electroosmotic flow which balances the flow perturbation generated by the DNA molecule. Hence, we can say that hydrodynamic interactions are screened and DNA is “free-draining” during free electrophoresis. Every segment feels the same friction from the fluid, so DNA molecules of different sizes will have the same mobility. We consider a DNA molecule subject to an external electric force E. The free electrophoretic mobility of DNA µ0 is defined as the ratio between the electrophoretic velocity v and electric field strength E. For a “free-draining” DNA bead-spring chain, the total friction coefficient of the chain is the sum of all the friction coefficients ζ of the beads. Likewise, the total electric charge of the chain is Q = Nb qb where qb is the effective charge per bead. Balancing the electric force Fe = QE = Nb qb E and the

19 friction force Fh = Nb ζv, results in the following center of mass velocity

v=

Nb qb E qb E = . Nb ζ ζ

(3.11)

Hence, the electrophoretic mobility of DNA in free-solution (µ0 ) is given by

µ0 ≡

qb v = . E ζ

(3.12)

This expression indicates that the electrophoretic mobility of DNA in free solution is independent of the length or the molecular weight of the chain. For a “free-draining” chain with diffusivity D, the friction coefficient is ζ = kB T /Nb D. Therefore, given the electrophoretic mobility, the effective charge per bead is

qb = ζµ0 =

kB T µ0 . Nb D

(3.13)

Using the diffusion coefficient of λ-DNA during free solution electrophoresis D = 0.55 µm2 /s, and the electrophoresis mobility obtained under the same condition µ0 = 4.2 ×10−4cm2 /Vs [102], qb is determined to be 178 electrons per bead, for the Nb = 11 bead-spring model used in this work.

3.5

Polymer-wall steric interaction

A potential needs to be devised to prevent penetration of beads through solid walls. However, if the potential is steep near a wall, such as that of a Lennard-Jones potential, a very small time step should be used. On the other hand, in the case of

20 soft potential, a bead can sometimes penetrate the wall. Jendrejack et al. devised a potential barrier for bead-wall interactions which is harder than the Gaussian soft potential used for bead-bead interaction and softer than typical Lennard-Jones potential [77], Uiwall =

   −2 3 Awall b−1 k δwall (hi − δwall ) , hi < δwall   0,

(3.14)

hi > δwall.

Here hi represents the perpendicular distance of bead i from the wall, δwall is the cut-off distance. The choice of the repulsive potential is arbitrary as long as it can prevent the bead from penetrating the wall with reasonably small time step. In Appendix A, we derive an effective repulsive potential between a polymer molecule and a hard wall using the ideal chain model.

21

Chapter 4 Bistability and field-driven dynamics of confined tethered DNA DNA is “free-draining” during electrophoresis as discussed in Chapter 3. As a consequence, it cannot be separated by free electrophoresis. A sieving medium, such as polymer gel or microfabricated structure, that induces a size-dependent mobility is necessary for DNA separations. Many microfluidic devices with well-defined structures have been proposed for DNA manipulation, especially for separation and stretching purposes using electrophoresis [45]. These devices are typically composed of periodic obstacle arrays. When an electrostatic field is introduced, a nonuniform field with steep gradient is generated in the device. Thus, understanding DNA dynamics in nonuniform electric fields is an important issue for the design of high performance mi-

22 crofluidic devices. The motivation of the work presented in this chapter1 is to build a general framework to analyze DNA dynamics in nonhomogeneous electrostatic fields in complex microfluidic geometries. A Brownian dynamics/Finite element method algorithm is presented in Appendix B. In the remainder of this chapter, we focus on its application to study the properties a novel class of soft nanomechanical control elements we proposed for microfluidic devices.

4.1

Introduction

A fundamental unit in many engineering systems is an element with two distinct states – on/off, open/closed, left/right etc. Of particular interest and utility in many cases is the property of bistability: each state is robust in the sense that a finite perturbation must be temporarily introduced to switch the system from one state to the other. Many classical macroscopic mechanical systems display bistability, and efforts have recently been made on many fronts to design, construct and analyze “soft” bistable systems on smaller and smaller scales, down to the molecular level. Steen and coworkers [70, 152] have constructed systems of pairs of droplets connected by a flow channel. This system exhibits bistability via the existence of situations where two configurations with minimal surface energy can arise, with one drop large and the other small, or vice versa. Switching between states can be achieved, for example, via pressure fluctuations or electroosmotic pumping of fluid between droplets. Groisman et al. [63] and Arratia et al. [10] have demonstrated microfluidic systems in which elastic liquids (solutions of long flexible polymer molecules) in complex geometries 1

This chapter is based on the publication: Y. Zhang, J. J. de Pablo, and M. D. Graham, Soft Matter, 5, 3694.

23 can exhibit bistability via symmetry-breaking flow instabilities. At the level of single flexible molecules in solution, bistability arises through coil-stretch hysteresis in sufficiently long DNA molecules in solution during extensional flow [126]. (We will use the term “bistable” loosely here, to denote situations where the energy barrier between two metastable states is greater that about kB T , where kb is Boltzmann’s constant and T is absolute temperature.) A similar phenomenon is predicted for a chain tethered on a no-slip wall, in the case where there is a stagnation point flow and the chain tether point is at the stagnation point [18]. Banavar et al. [15] used simulations of model chain molecules to predict bistability in the conformations of the molecules. In particular, they predict for a certain case the coexistence of a single helix and a dual helix in which the chain folds in half and self-interacts to form a double helix. Our interest here is in bistability of polymer chain conformations associated with the combination of tethering and confinement within a micro- or nanoscale device. Partially confined configurations of polymer chains, where different parts of a chain experience different degrees of confinement, give rise to interesting phenomena of entropic origin. The entropy difference between a portion of polymer chain in confined space and in open space gives rise to an entropically induced recoil force, which tries to pull the whole chain into the open space. Using a nanochannel/microchannel device, Mannion et al. estimated this entropically induced force of double-stranded (ds-) DNA to be about 102 − 103 f N[92]. According to the force versus extension curve of ds-DNA, a force of this magnitude will generate a fractional extension of the molecule of about 0.5 [94]. In a related calculation, Bickel et al. determined the force required to tether an end of an ideal chain onto a hard wall, which is also a force due to the decrease of available number of configurations [20]. This tethering force converges to

24 a constant value kT /lk as chain length diverges. For DNA with a Kuhn length of lk ≈ 100nm, this force is about 40 f N at room temperature. Partial confinement has been exploited to manipulate a DNA molecule using a nanoslit device (50 - 100 nm across), on which one wall had square nanopits (100 nm deep and 100 - 500 nm wide)[116]. Because the confinement is weaker in the regions of the slit where there is a pit, DNA segments are preferentially found there. This entropic well tends to pin sections of the chain in these regions. When multiple pits are present (and not too far apart), stretched sections of DNA will span the more confined regions between pairs of pits. Binder and coworkers have used Monte Carlo simulations and theory to study a tethered polymer in good solvent compressed by a circular disk centered a distance H above the tether point [95, 71]. For moderate compression, the polymer is “imprisoned” underneath the disk adopting a quasi-2D random walk configuration by forming blobs with diameter equivalent to the disk height H. When a certain height Himp , part of the chain “escapes” this imprisonment and forms a “stem-and-flower” configuration where the blobs underneath the disk are stretched into an extended configuration by the entropic force exerted by the escaped parts. Using Monte Carlo simulations, they found that there is a regime of heights Hesc ≤ H ≤ Himp where the chain exhibits “two-phase coexistence” (bistability) of the escaped and the imprisoned configurations. Both simulation results and analytical theory indicate that the escape transition is a first-order phase transition. This provides an example of bistability created by a combination of tethering and partial confinement, specifically by the competition between the entropy gain of the escaped segments and the entropy loss and energy penalty of extending the portion of the chain that remains trapped under the disk.

25 Partially confined configurations also lead to interesting dynamics [103, 99]. Trapped within a two-dimensional array of spherical cavities interconnected by circular holes, short linear DNA strongly localize in cavities and only sporadically “jump” through holes [103]. To jump out of a cavity, small DNA segments penetrate into the neighboring cavity first and form a partial confined configuration. The fluctuation of the probing segments beyond a threshold leads to an abrupt jump of the entire molecule. This is in accordance with previous theoretical studies of entropic barriers by Muthukumar [100, 99]. In the present work, we examine, using Brownian dynamics simulations of a coarsegrained model for long flexible DNA molecules, the possibility of generating bistable behavior using long flexible linear polymer molecules end-tethered in a confined geometry. Fig. 4.1(a) shows a very simple example of what we will call an “entropically bistable” system. (The system studied by Binder and coworkers [95] is also in this class.) Here a polymer molecule is tethered within a small pore of width 2d connecting two open regions. For a sufficiently long chain (Rg ≫ d, where Rg is the radius of gyration of the chain in the absence of confinement), one clearly expects that the most entropically favorable situations will be those where the entire chain is in one open region or the other. Thermal fluctuations or transient application of an electric field (if the chain is a polyelectrolyte) or flow field might drive the chain from one region to the other over an entropic barrier. One realization of this type of system would be the fully three-dimensional situation where the pore is a hole connecting two half-spaces, as shown in Fig. 4.1(d). A more experimentally tractable version, fabricated using nanolithography techniques, is the “quasi-two-dimensional” (Q2D) situation shown in Fig. 4.1(c), where top and bottom bounding walls exist. Both

(b)

y

(a) x

y

z thermal noise

x

z

2d or

2d

transient field

no field

l

high field

medium field

(d)

(c)

r

y

h=2d

x z

y

x z

Figure 4.1: Schematic representations of soft nanomechanical bistable elements. (a) Pore-crossing geometry: an entropically bistable system. (b) Pore-entry geometry: a competitively bistable system. (c) Quasi-2D porecrossing geometry. (d) 3D pore-crossing geometry.

26

27 situations will be studied here. We also consider the possibility of systems where the origin of bistability is somewhat more subtle. Consider the system shown in Fig. 4.1(b), comprising a polyelectrolyte end-tethered at the mouth of a long pore. In the absence of an electric field, the most likely conformations of the chain lie outside the pore (left). If a sufficiently large field is applied pulling the chain into the pore, entropy will be overcome by electrostatic energy and the chain will reside entirely in the pore (right). We hypothesize that at intermediate field strengths, a nontrivial competition between entropy and electrostatic energy will allow the possibility of two stable states (center), one where the chain is mostly outside the pore, in a high entropy state, and the other where the chain is mostly inside the pore, in a low entropy state. This system will be called “competitively bistable”. We will study the quasi-two-dimensional version of this geometry. Many fundamental issues arise when considering these systems. What is the nature of the transition states in these systems? How high are the free energy barriers between states? How do these depend on the length of the chain and, in the Q2D case, the degree of lateral confinement? What are the dynamics of transitions between the states? How will those dynamics differ in the electrostatically driven and flow driven cases? To what extent do the dynamics of these systems follow Kramers-type kinetics? This initial report will only begin to touch on some of these issues. Potential applications of these systems can also be considered. A chain tethered in a micro/nanochannel under sufficiently confining conditions would behave somewhat like a porous medium, resisting fluid flow and blocking the transport of sufficiently large solutes. This property in combination with entropic bistability might be ex-

28

(a)

(b) z Transient Field

x

z y

x

y

Transient Field

Figure 4.2: Possible applications of tethered polymer molecules as control elements in micro/nanofluidic devices. (a) Switch. (b) Gate.

ploited as shown in Fig. 4.2. In the left image, a polyelectrolyte chain is configured so that it blocks transport of the round solute particles in the left channel, while allowing transport of the square ones, or vice versa. Switching would occur through transient application of an electric field in the z direction. In the right image, an entropically bistable geometry serves as a gate that can be opened or closed by a transient electric field.

4.2

Polymer model and simulation approach

We use Brownian dynamics simulation of a model of long ds-DNA molecules to study the statics and dynamics of a tethered chains in confined systems. Following previous work [76, 74], a double-stranded DNA molecule is described by a bead-spring chain model composed of Nb beads of hydrodynamic radius a = 77 nm connected by Ns = Nb − 1 entropic springs. Each bead represents a DNA segment of 4850 base pairs, i.e., Nb = 11 corresponds to a stained λ-DNA, which has a contour length of 22 µm

29 (a)

(b)

Pore crossing

Pore entry

y

l=2 a

x

z

ze

2d=14 a

ze

2d=14 a

Figure 4.3: Top view of the simulation domains. Black areas are impenetrable walls; the free end of the chain is a filled circle. In quasi-2D cases, the tethered bead-spring chains are bounded in x direction by walls at x = −d and x = +d. (a) Pore-crossing geometry. (b) Pore-entry geometry.

and radius of gyration of 730 nm. An effective charge of Q = 178 e is assigned to each bead, which is calculated based on the free solution mobility of DNA [142, 102]. The length unit in this work is a = 77 nm. The springs connecting the beads obey a worm-like chain force law [94] Fsij

  |rj − ri | −2 4|rj − ri | rj − ri kB T (1 − ) −1+ ] . = 2bk Nk,s bk Nk,s bk |rj − ri |

(4.1)

Here, bk is the Kuhn length for DNA and Nk,s = 20 is the number of Kuhn length per spring. The physical confinement is taken into account through an empirical bead-wall repulsive potential of the form

Uwall = i

   −2 3 Awall b−1 k δwall (hi − δwall ) , hi < δwall   0

hi > δwall .

(4.2)

30 Here hi represents the perpendicular distance of bead i from the wall, δwall is the 1/2

cut-off distance. In this work, we choose Awall = 50kB T /3 and δwall = bk Nk,s /2 = 3.01a = 0.24 µm. In this work, we ignore hydrodynamic interactions (HI) between chain segments and assume that the chain is ideal. The force balance on the beads leads to a stochastic differential equation[105]: 1 dR = Fdt + ζ

s

2

kB T dW, ζ

(4.3)

where R is the vector containing bead positions ri , ζ is the friction coefficient of the bead and F is the vector of non-hydrodynamic and non-Brownian forces. The components of dW are obtained from a real-valued Gaussian distribution with mean zero and variance dt. We use a standard (stochastic) Euler scheme[105] for timeintegration. In the pore crossing system (Fig. 4.3 (a)), we study the dependence of the free energy barrier on chain length and confinement in both a quasi-2D and 3D geometries. The wall thickness is 2 a and the pore size is 14 a. For the quasi-2D case (Fig. 4.1 (c)), the chain is confined in the x direction, and the distance between the two impenetrable planes are h = 14 a. For the 3D case (Fig. 4.1 (d)), the chain is unbounded in both x and y directions. In the pore entry system (Fig. 4.3 (b)), we consider only the quasi-2D case. The chain is tethered to the center of the entrance to a square channel which runs along the +z axis. An electrostatic field is applied along the −z direction. In the present model electrostatic interactions between beads are neglected, consistent with behavior in a high ionic strength solvent, where electrostatic interactions are screened. Elec-

31 troosmosis of counterions screens hydrodynamic interactions, as does confinement, justifying their neglect in the present simplified model. The electric field strength E at z = −∞ is E = E ∞ ez . The relative field strength is denoted as ER = αE ∞ /E0 , where E0 = kB T /aQ = 18.74 V /cm is the field strength unit and α is a constant to make ER order O(1). When ER = 1, the corresponding E ∞ = 1.23 V /m, and the electric force f E = E ∞ Q = 3.51 × 10−2 fN. Even inside the channel where the electric field gradient is around ten times larger than E ∞ , the electric force is very small compared with thermal agitation kB T /a = 53.4 fN. The electric field is obtained by solving the governing Laplace equation with no flux boundary conditions on the channel walls and Dirichlet boundary conditions at the left and the right edges of the domain. The commercial PDE solver COMSOL, based on the finite element method (FEM), is used. When d = 7a, the simulation boxes are y × z = 200a × 100a and 14a × 50a for the outer chamber and the channel, respectively. As the field gradient becomes uniform and constant going along the channel, we can use a short simulation box for the channel. We use the Lagrange-quadratic square elements (1a × 1a) provided by COMSOL. As the polymer model is coarse-grained to a length scale of a, this resolution of the field is adequate.

4.3 4.3.1

Results and discussion Pore-crossing geometry

We consider first the pore-crossing geometry of Fig. 4.3(a). Only the equilibrium behavior in this system will be studied. Fig. 4.4(a) is a typical time series plot of the z component of the free end ze of a chain with Nb = 15 in the 3D case. This plot

32

(a)

4

4

(b)

( m)

2

0

z

z

e

0

e

( m)

2

-2

-4 0

t

-2

w

500

1000

1500

-4 1070

2000

1080

t (s)

1090

1100

t (s)

(c)

0.2

w

(t )

0.3

t

/ exp(-t /

w

w

0.1

0.0 0

1000

2000

3000

4000

5000

t (s) w

Figure 4.4: Simulation results for an ideal chain tethered to the center of the pore in a 3D pore-crossing geometry (Fig. 4.3(a)). (a) Typical time series plot of the z component of the free end of a tethered chain with Nb = 15 and d = 7a. (b) Blowup of a segment of the time series plotted in (a). (c) Waiting time distribution and best fit to an exponential distribution.

33

(a)

2

(c)

2

1

1

1

0

0

0

−1

−1

−1

0

1

2

−2 −3 −2 −1

3

2

(e)

0

1

2

−2 −3 −2 −1

3

2

(f)

1

1

0

0

0

−1

−1

−1

0

1

2

3

−2 −3 −2 −1

0

1

2

3

1

2

3

0

1

2

3

2

1

−2 −3 −2 −1

0

−2 −3 −2 −1

y( mm )

−2 −3 −2 −1

(d)

(b)

2

z( m m) Figure 4.5: Snapshots of a crossing event in the quasi-2D pore-crossing geometry. Time interval between frames is 0.5s.

clearly shows that ze fluctuates around two metastable states. Fig. 4.4(b) is a blowup of the time interval 1070-1100 s. The chain end stays on one side of the wall for a time interval tw , which we will call the waiting time, and then crosses to the other side. Fig. 4.5 shows the dynamics of a typical crossing event. The free end of the tethered chain first diffuses to the pore and then through to the other side of the wall. We notice that this step of end bead crossing does not necessarily cause a crossing of the whole chain. A successful crossing event occurs when the number of beads on one side of the wall changes from Nb − 1 to 0 or vice versa. The time points tic when a crossing completes are recorded. We define the waiting time as the time between two adjacent crossing events tiw = ti+1 − tic . As seen in Fig. 4.4(a), the crossing time c is very small compared to the waiting time, which is well-fit to an exponential with chain-length-dependent decay rate 1/τ , as shown in Fig. 4.4(c). This result indicates

34 that, to a good approximation, the crossing events are independent of one another – crossing is a Poisson process. To gain further information about the crossing process and the transition state, we consider the joint behavior of the z-component of the chain end, ze and the zcomponent of the center bead, zm (where m = Nb /2 and (Nb + 1)/2 when Nb is even and odd respectively). The joint probability density for these two variables is denoted ρ(ze , zm ). Fig. 4.6(a) shows − ln ρ(ze , zm ) for a chain with Nb = 13 chain in a 3D pore-crossing geometry. It clearly demonstrates that there are two metastable states (maxima in ρ(ze , zm )), in which both ze and zm are on the same side of the wall, and two transition states when the end bead and the center bead are located on different sides of the wall. The transition states are in accordance with what we observed in the snapshots shown in Fig. 4.5. A simpler representation of the free energy surface explored by this system is given simply by the probability density for the z-component of the end bead ρ(ze ). Unlike the joint distribution ρ(ze , zm ), this representation is too low-dimensional to provide insight into the structure of the transition state. On the other hand, it is simple to construct and visualize, and a simple potential of mean force (PMF) can be constructed as  ρ(ze ) , F (ze ) = F − kB T ln ρ∗ ∗



(4.4)

where F ∗ and ρ∗ are arbitrary reference values. In Fig. 4.6(b) we show ρ(ze ) (dashed circle) and corresponding F (ze ) (solid line) when F ∗ = Fwell = 0 (Nb = 20) in a quasi-2D system . As expected, the PMF is composed of two energy wells which are separated by a energy barrier – the difference in energy between these is denoted ∆Fmax . Fig. 4.6(c) shows a plot of the energy barrier ∆Fmax /kT as a function of

35

2 -1.5

z ( m)

1

-3.5 -3.0

m

-4.0

0

-2.5 -3.5

-1 -4.0

-2

-2.0

3

(b)

2

e

-5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.50 0

F(z )/kT

3

(a)

1

F

-0.50

max

-3.0

-1.0

-3 -4

-2

0

2

0 -3

4

-2

-1

7

F

max

/kT

5

(d)

4

F

3 2 1 0

3D

2

3

1000

a=2.58 0.12

4.16 0.14

max

kT

a ln(N ) b

(s)

6

1

e

e

(c)

0 z ( m)

z ( m)

100 b

3D

a=1.41 0.03

Q2D

10

Q2D

2.45 0.05

10

15 N

b

20

25 30

10

15 N

20

25

30

b

Figure 4.6: Simulation results for an ideal chain tethered to the center of the pore in quasi-2D and 3D pore crossing geometries (Fig. 4.3(a)). (a) The negative of logarithm of the joint probability density function of the z-component of the end bead ze and the z-component of the center bead zm for a Nb = 13 chain in 3D. (b) Reduced potential of mean force ∆F (ze )/kB T for a Nb = 20 chain in Q2D. (c) Semilog plot of reduced energy barrier vs. chain length when d = 7a (d) Log-log plot of mean first passage time τ v.s. chain length when d = 7a. The error bars are about the same size as the symbols.

36 Nb for the 3D and quasi-2D cases. For the range of Nb considered it appears that ∆Fmax /kT = a ln Nb . Fig. 4.6(c) shows the dependence of the mean waiting time τ on molecular weight. Over the limited range of molecular weights studied, the results appear to follow a power law: τ = τm Nbα . If the transition process follows Kramers kinetics [160], one would expect that τ = τ0 exp(∆Fmax /kb T ), which in the present case would imply that τ0 = τm Nbα−a . If hydrodynamic interactions were included in the model, one might expect α to decrease slightly, based on simulation results for translocation of a free polymer chain through a pore ([8, 73, 56, 54]). Since a is determined from the free energy landscape of the system, it is unaffected by the presence or absence of hydrodynamic interactions in the model. On the other hand, the pore geometry (size, shape) might affect both these parameters; in particular changing the pore size will change the free energy landscape and thus the exponent a. The nature and origin of the observed exponents will be the subject of future work; we note that extension of Kramers’ theory of barrier crossing of a point particle to a polymer system is a topic of some current research [106, 129, 128]. Finally, we observe that the simulation results indicate that confinement from 3D to quasi-2D greatly reduces the energy barrier and leads to a larger transition rate (Fig. 4.6(c) and (d)). As the crossing rate is inversely proportional to the free energy barrier, which is closely related to the number of accessible chain configurations, this is not surprising by considering that the number of accessible chain configurations in the quasi-2D case is greatly reduced compared to that in the 3D case.

37

4.3.2

Pore-entry geometry

The previous section reported a system where geometry alone determines the free energy landscape of the tethered polymer. The present section turns to a case where entropy and electrostatic energy compete, and address the possibility that bistability of a tethered chain can arise from this competition. The situation under consideration is the pore-entry geometry of Fig. 4.3(b). Only the quasi-2D case will be considered. In the absence of a field, a chain with Rg ≥ d would not be expected to sample the pore, while if the field is sufficiently large it is expected that the entire chain will reside in the pore. The regime of interest is chain lengths such that Rg ≥ d and intermediate values of the field strength. Fig. 4.7(a) shows a time series of ze for the case Nb = 10 and ER = 5. As in the pore-crossing geometry, bistability is observed, which can be seen more clearly from the joint probability density function ρ(ze , zm ) illustrated in Fig. 4.7(b). The path connecting the two metastable states indicates, again, the transition state. Fig. 4.8(a) shows PMF curves based on ρ(ze ) as functions of ER when Nb = 10. As the electric force is small in the outer chamber, it does not significantly contribute to the energy of the chain when ER . 10.0 and the energy profile outside the channel is insensitive to the change of ER . The energy well inside the channel is created by a balance between the entropic force and the electric force. We denote the critical field strength above which this energy well appears as Elow (≈ 2.2 for Nb = 10). As ER increases further, the free energy barrier for entry to the pore weakens, and eventually vanishes once ER exceeds a field strength denoted Ehigh When Elow < ER < Ehigh , the system shows competitive bistability that is created by the interplay between entropy and electrostatic energy. Fig. 4.8(b) shows the positions of the minima in Fig. 4.8(a) as

38

4

(a)

3

e

z ( m)

2 1 0 -1 -2 0

20

40

60

80

100

t(s)

(b) 4

-4.00 -3.50 -3.00 -2.50 -2.00 -1.50 -1.00 -0.500 0

3

1

m

z ( m)

2

0

-1 -2 -3 -3

-2

-1

0

1

2

3

4

z ( m) e

Figure 4.7: Simulation results for a tethered ideal chain in the quasi-2D pore-entry geometry (Fig. 4.3(b)) for the case Nb = 10, ER = 5, d = 7a. (a) A typical time series plot of the z component of the free end of a tethered chain. (b) The negative logarithm of the joint probability density function ρ(ze , zm ).

4

(a)

(c)

3 2 1

e

F(z )/kT

0

E1 E2.2

-2

E3

-3

E4 E5

-4

z

E7

-5

(e)

e0

1

0

0

−1

−1

−2 −2 2

−1

0

1

2

-1

0

1

2

3

4

3 OUT

h

IN

2

(g)

−2 −2

3

(f)

1

0

0

−1

−1 −1

0

1

2

3

(h)

0

1

2

3

−1

0

1

2

3

−1

0

1

2

3

2

1

−2 −2 2

−1

−2 −2 2

y( m m)

-2

e

f

1

d

e0

( m)

2

1

E9

z ( m)

(b)

(d)

E0

-1

-6 -3

2

E

z

E

high

low

0

-1

-2

0

2

g

e

c

4

6 E

R

8

10

12

14

1

1

0

0

−1

−1

−2 −2

−1

0

1

2

−2 −2

3

z( m m)

39

Figure 4.8: Simulation results of an ideal chain tethered to the center of the entrance in a quasi-2D pore entry geometry, with Nb = 10. (a) Potential of mean force as function of last bead position ze and reduced electric field gradient ER . Top to bottom curves correspond to ER = 0, 1.0, 2.2, 3.0, 4.0, 5.0, 7.0, and 9.0, respectively. (b) Most probable end bead positions ze0 (i.e., maxima of ρ(ze ) or minima of ∆F (ze )) as a function of ER for state ”IN” and state ”OUT”. The error bars are almost the same size as the symbols. (c)(d),(e)(f),(g)(h) are snapshots at ER = 2.3, 5, and 9, respectively.

40

(b)

(a) 14

1.0 0.9

12

E

0.8

low

E

0.7

high

P

in

8

E

R

10

State "IN"

6

N =10 b

0.5

N =15 b

0.4

N =20

0.3

4 State "IN"+"OUT"

0.04

0.06 1/N

0.08

0.10

b

0.1

State "OUT"

0.02

b

N =25

0.2

2 0 0.00

0.6

0.12

0.0 0.0

0.5

b

1.0

1.5

2.0

2.5

3.0

E /E R

1/2

Figure 4.9: (a) Phase diagram showing the region in parameter space where the chain shows bistability in the pore-entry geometry. (b) The probability to find the chain inside the channel Pin as function of ER normalized by E1/2 . functions of ER – the regime of bistability is clearly seen. Fig. 4.8(c) shows typical conformations of the chain at various points on Fig. 4.8(b). Finally, we address the phase transition behavior of the tethered polymer molecule in the competitively bistable system. In Fig. 4.9(a), the shaded area shows the parameter regime in ER and 1/NB where bistability occurs. Fig. 4.9(b) shows the probability R∞ Pin = 0 dze ρ(ze ) of the chain end residing inside of the channel, as a function of ER

for various Nb . We have normalized ER with E1/2 , the (molecular-weight dependent)

field strength where Pin = 0.5. From these plots it appears that as N → ∞, the slope diverges and the region of bistability vanishes, indicating a first-order phase transition. Similar behavior was noted in the escape transition system studied by Binder and coworkers [95].

41

4.4

Conclusions

Simulations predict that flexible polymers end-tethered in a pore within a confined geometry can display multiple free energy minima separated by substantial barriers – bistability. For a tether point in a pore between two large chambers, bistability is expected based on simple considerations of entropy and geometry -we denote this as “entropic bistability”. In this situation dynamic simulations suggest that the crossing of the chain end through the pore is the dominant mode of transition between the two states, and consideration of the joint probability density function of the position of the chain end and the chain midpoint indicates the presence of a transition state where the chain end is on one side of the pore, while the chain midpoint is on the other. For a chain end-tethered at the mouth of a long pore and subjects to a field that tends to draw the chain into the pore, we predict the existence of a second kind of bistability, arising from the interplay of conformational entropy and electrostatic energy; we call this “competitive bistability”. In particular, simulations show the presence of a regime of field strength in which there are two free-energy minima, one corresponding to a chain that is largely outside the pore in an entropically favorable state, and one where the chain is largely inside the pore, in an electrostatically favorable state. Both fundamental and technological question are raised by these results. The systems under consideration here comprise model systems for activated rate processes in systems with many degrees of freedom and provide an opportunity to compare experiment, theory and simulations for this class of systems. Of particular interest may be the difference in dynamics when transition from one state to another is driven by fluid flow rather than an electrostatic field. These systems or variations of them might also have applications as electrically or fluidically actuated nanomechanical or

42 nanoelectromechanical elements, such as valves or switches.

43

Chapter 5 Tethered DNA dynamics in shear flow In Chapter 5 and Chapter 6, we discuss flow-driven DNA dynamics. To study the effects of solid impenetrable walls on the dynamics of a nearby DNA molecule, we start by re-examining a problem with simple geometry: the cyclic dynamics of a tethered DNA molecule in shear flow where the chain is grafted to a flat wall. We compare three simulation methods: Brownian Dynamics (BD), the Lattice Boltzmann Method (LBM), and a recent Stochastic Event-Driven Molecular Dynamics (SEDMD) algorithm.1 We focus on the dynamics of the free end (last bead) of the tethered chain and we examine the cross-correlation function (CCF) and power spectral density (PSD) of the chain extensions in the flow and gradient directions as a function of chain length N and dimensionless shear rate Wi. Extensive simulation results suggest 1

This chapter is based on the publication: Y. Zhang, A. Donev, T. Weisgraber, B. J. Alder, M. D. Graham, and J. J. de Pablo, J. Chem. Phys. 130, 234902. LBM and SEDMD calculations were performed by our collaborators at the Lawrence Livermore National Lab.

44 a classical fluctuation-dissipation stochastic process and question the existence of periodicity of the cyclic dynamics, as previously claimed. We support our numerical findings with a simple analytical calculation for a harmonic dimer in shear flow given in Appendix C.

5.1

Introduction

The interaction of polymer molecules with fluid flow has been studied both theoretically [130, 21] and experimentally [109, 108, 136, 135, 57] for several decades. The behavior of polymer chains in flow is determined by an intricate interplay between the flow gradients, chain elasticity, thermal fluctuations, and the physical confinement [78, 75]. The dynamics of tethered polymer molecules (“polymer brushes”) in shear flow has received considerable attention due to its relevance to diverse important applications, such as colloidal stabilization, surface adhesion, and lubrication [96]. In contrast to previous work on the collective motion of polymer brushes [96, 55], Doyle et al. studied the dynamics of a single tethered DNA in uniform shear flow using fluorescence videomicroscopy [46]. Enhanced temporal fluctuations in the chain extension were observed, and were attributed to the coupling of advection in the flow direction and diffusion in the gradient direction. A cyclic dynamics mechanism (Fig.5.1), closely related to the tumbling dynamics of a free polymer molecule in shear flow [115, 57, 127], was proposed based on results from Brownian dynamics simulations. No peaks were observed in the calculated power spectral density (PSD) of the DNA extension in the flow direction, and the authors therefore suggested that the cyclic dynamics of a tethered chain in shear flow is aperiodic. An important

45 physical question is whether there is a characteristic timescale associated with the cyclic motion that is distinct from the internal relaxation time of the chain. Several computational studies revisited the problem of a tethered chain in shear flow by looking at different variables relating to both the flow and gradient directions, such as extensions along both flow and gradient directions [37], polymer orientation angle defined through the gyration tensor [127], and angle between the wall and the vector joining the tethering point to the center-of-mass of the chain [62]. The cross-correlation functions for such variables exhibit signatures of the proposed cyclic motion in the form of peaks at non-zero delay time. Because of the particular choice of variables in Ref. [62], the lack of such peaks at small Weissenberg numbers was attributed to the existence of a critical Weissenberg number; as demonstrated in Ref. [42], choosing a different sets of variables shows that the signature peaks exist even at small shear rates. In Ref. [62], the position of the peak in the studied cross-correlation functions was interpreted as a characteristic cycling time, and it was found to be a fraction of the relaxation time of the polymer chain. In Ref. [127] the tumbling motion of a free polymer chain in shear flow was studied experimentally and computationally, and wide peaks were found in the Fourier spectra of the time series of the angle of the chain relative to the flow direction. These peaks were identified as evidence of periodic motion of the tumbling molecule. The characteristic tumbling time (period) was extracted from the position of the peak in the spectrum and was found to be in good agreement with the experimentally-measured tumbling frequency. These finding for a free chain in shear flow inspired similar studies of a tethered chain, and similar observations of periodic motion with a characteristic period about an order of magnitude larger than the relaxation time were reported [127, 37, 38]. In

46

x

y

z

1)

2)

3)

11 1

2

4)

Figure 5.1: Snapshots taken from a simulation run with Wi = 5 to show tethered DNA dynamics. The beads are labeled from 1 to 11 as shown. Cyclic motion mechanism proposed by Doyle et al. is composed of four stages: 1) (Re)coiling; 2) initiating; 3) stretching; 4) rotating [46] several later studies of the tumbling emotion of a free polymer chain in shear flow, experimental results [57], numerical simulation [115], and theory [34, 156] all suggest that the intervals between successive tumbling events are exponentially-distributed with a decay constant equal to the relaxation time of the chain. Such an exponential tail implies that the tumbling events occur as a Poisson-like process, which is aperiodic despite the existence of a characteristic timescale (frequency of repetition). If the tumbling events of a free polymer in shear flow is aperiodic, intuitively, adding of a wall to break the symmetry of the motion should leave the dynamics of a tethered chain aperiodic. Different authors use the terms “cyclic” (repetitive) and “periodic” with differ-

47 ent meanings, and it is therefore important to give our definitions. Periodicity is the quality of occurring in regular time intervals (periods). Periodic motion has correlation functions that are (possibly damped) oscillatory functions, and spectra that have sharp peaks. Noise (fluctuations) and the associated dissipation will always broaden any peaks that are related to underlying deterministic periodic motion (and consequently, exponentially dampen the oscillations in the real-space correlation functions). As an example, for a rigid spheroid in shear flow, there is indeed periodic motion (Jeffery’s orbits) in pure shear flow. Adding fluctuations, when they are small, is expected to preserve that but introduce some broadening of the spectral peaks [88]. In contrast, a cycle usually means a process that eventually returns to its beginning and then repeats itself in the same sequence. The end-to-end tumbling of a single polymer molecule in shear flow provides a relevant example. In this paper we analytically calculate the power spectrum for a tethered dimer in shear flow, and find an exponentially-decaying cross-correlation function that has the relaxation time as the only characteristic timescale. More importantly, this analytical example shows that the power spectrum can exhibit a wide peak at small frequencies without any underlying periodic motion, and that the location of a maximum in the PSD is not necessarily an indication of a new timescale. The analytical results for a tethered dimer are consistent with our numerical observations for longer tethered chains in shear flow. Therefore, our investigations do not confirm the existence of periodic motion with a period distinct from the relaxation time of the chain, as previously suggested in the literature [127, 37, 38]. We apply three different solvent representations to the same problem of a tethered chain in shear flow. The models are an implicit solvent (Brownian Dynamics

48 [74], BD), a continuum solvent (Lattice-Boltzmann [5], LB), and a particle solvent (Direct Simulation Monte Carlo [42], DSMC). Such comparison between widely differing methods on the same problem is important as a validation of their range of applicability. It is also important to compare the computational performance of the different methods. In this work, different chain representations and boundary conditions make a direct quantitative comparison impossible. Specifically, the BD polymer is a worm-like chain representative of semi-flexible DNA, in the LB simulations it is a flexible chain of repulsive spheres, and in the DSMC solvent it is a flexible chain of hard spheres. However, we can access the importance of the details of the chain model, and in particular, of chain elasticity, and thus test the widely used assumption that the dynamics scales with the Weissenberg number Wi = γτ ˙ , where γ˙ is the shear rate and τ the chain relaxation time, independent of the details of the model. For this particular problem of a tethered chain in shear flow, we find good agreement between the different methods. A general discussion of the wide range of techniques for modeling the hydrodynamics of polymer chains in solution is given in Section 5.2. Further details about the three specific techniques we use in this paper to study the tethered polymer problem are given in Section 5.3. In Section 5.4 we present our results, and finally, in Section 5.5 we give some concluding remarks.

49

5.2

Discussion of Methods for Hydrodynamics of Polymer Solutions

In this Section we give a brief overview of various methodologies for modeling hydrodynamics of soft matter systems, notably, polymer solutions (see also review by Duenweg and Ladd [47]). The various methods for computational hydrodynamics of polymer solutions can be divided in two major categories. The first are purely continuum methods that use constitutive equations for the polymer solution. These models only apply at macroscopic scales, when the number of polymer chains in an elementary fluid flow volume is large, so that statistical averages of the chain conformations can be used as parameters in constitutive models of the time-dependent stress as a function of the strain rate history. The construction of such constitutive models is ad hoc and rather difficult in situations where conformations of the chains couple to an unsteady flow, as, for example, in the problem of turbulent drag reduction. Additionally, such continuum methods do not apply to situations where the dynamics of individual chains are of interest, such as a DNA molecule flowing through a micro-channel or DNA translocation through a pore. The second major category of methods explicitly simulates the motion of each polymer chain using some form of molecular dynamics. The simplest chain model is a dumbbell. Multi-bead representations of the chains are capable of complex chain conformations but require models for the bead-bead interactions. Such details of the polymer model are important for both physical fidelity and computational efficiency. For example, preventing chain-chain crossing can require stiff interactions for excluded volume terms, which in turn can lead to small time steps. There are also two major types of algorithms

50 for dealing with the solvent. One represents the solvent implicitly, and the others use an explicit solvent. An implicit solvent is most efficient computationally, however, it can only be used if the fluid flow in the absence of the polymer is known analytically or can easily be pre-computed numerically (e.g., stationary flow), and if the polymer chains themselves do not alter the background flow.

5.2.1

Implicit Solvent: Brownian Dynamics

The most widely used implicit-solvent algorithm is Brownian dynamics [74], described in more detail in Section 5.3.1. The method involves solving first-order differential equations of motion for the positions of the beads with additional forces due to the presence of the solvent. These solvent forces can be separated into a deterministic portion, for which a (linear) analytical approximation is used, and a stochastic portion, which is assumed to be white noise. The fluctuation-dissipation theorem is used to set the magnitude of the stochastic forcing. Brownian dynamics relies on several assumptions usually valid in microfluidic applications. The first assumption is that of small Reynolds number laminar (usually stationary) flow adequately described by a linearized Navier-Stokes equation. The second assumption is that hydrodynamic fields develop infinitely quickly relative to the rate at which the polymer conformation changes, so that a quasi-stationary approximation can be used to describe the perturbation of the flow field induced by the motion of the beads. This approximation leads to Stokes friction on single beads, as well as hydrodynamic interaction pairwise terms approximated with a long-range Oseen tensor as derived through an asymptotic (t → ∞) analysis for point particles. The free-draining approximation of Brownian dynamics neglects these pairwise hydrodynamic interactions. The in-

51 clusion of pairwise hydrodynamic interactions leads to a matrix formulation of the fluctuation-dissipation theorem and therefore factorization of a matrix of the size of the number of beads is required at every time-step. Various numerical tools have been devised to avoid dense factorization [76, 132, 74, 67], thereby reducing the cost of a single time step in Brownian dynamics with hydrodynamic interactions. Brownian Dynamics should be distinguished from Langevin dynamics, in which second-order (Newton’s) equations of motion are used for the beads, that is, both the bead velocities and positions are included as explicit degrees of freedom (but the solvent is still implicit) [122]. This assumes that there is a large separation of timescales between the fluid degrees of freedom and the velocities of the beads, which is in fact only true if the beads are much denser than the solvent. Furthermore, a much smaller timestep necessary to resolve the faster dynamics (relaxation) of the bead velocities. Therefore, Langevin dynamics finds its use only when the solvent is represented explicitly, so that calculating the friction and stochastic forces no longer requires factorization of the mobility tensor. An important advantage of Brownian dynamics is that it simulates the limit of zero Reynolds number exactly. It can also often exactly account for simple boundary conditions (e.g., flow in an infinite half plane) without resorting to approximations that truncate the flow field to a finite domain, such as the commonly-used periodic boundary conditions. Brownian dynamics is relatively easy to implement, however, complex boundary conditions, such as indentations or bumps on walls, requires care so that analytical approximations to the Oseen tensor that preserve the positivedefiniteness of the diffusion tensor [69]. While the computational cost can rise rapidly as the number of beads is increased when direct implementations are used, novel

52 schemes can be used to truncate the long-range hydrodynamic interactions and yield a linear dependence on system size, similarly to the handling of electrostatic interactions in spectral [132] and multipole methods [67].

5.2.2

Explicit Solvent: Continuum Methods

In order to capture the bi-directional coupling between the motion of the polymer and the flow around it, it is necessary to explicitly represent the solvent. The first level of approximation is to use a continuum description of the solvent assuming the applicability of the Navier-Stokes (NS) PDEs at small length scales. Typically an incompressible assumption is made, which is appropriate at sufficiently low Mach numbers if acoustic waves are not of interest. Additional approximations such as linearization or an iso-thermal approximation may be appropriate. The time-dependent (unsteady) NS equations can be solved by any of the numerous existing CFD algorithms, including explicit, implicit, or semi-implicit algorithms of varying level of complexity [131, 146, 12]. One of the advantages of the PDE formulation over particle methods is the ability to use powerful adaptive mesh resolution techniques that allow coarsening of the mesh away from the region of interest, here polymer chains. However, the case of complex boundary conditions such as needed, for example, in the handling of moving beads or flow through porous media, presents difficulties. An alternative to solving the Navier-Stokes PDEs is to use the Lattice-Boltzmann (LB) method [150], as discussed in Section 5.3.2. It requires small time steps limited by CFL-type conditions, however, each of the time steps is efficient. Recently, so-called entropic LB schemes have been developed that posses a discrete H-function, resulting in unconditional numerical stability even at high Reynolds numbers [23]. LB has been

53 found competitive with NS solvers in many situations and has the further advantage that it is based on kinetic theory and allows a more detailed level of description than NS. An important advantage of LB solvers is also their ability to handle complex boundary conditions. Recently, Chen et al. have provided a detailed comparison between BD and LB simulations on a DNA model that shows that the LB method provides a reasonable description of the results of more precise BD simulations at low Reynolds numbers [32]. Thermal Fluctuations Most continuum fluid dynamics methods are deterministic and thus do not include internal fluctuations of the hydrodynamic fields. Fluctuations become more important the smaller the length scale of interest, and are crucial for polymer flows. Including thermal fluctuations in a continuum formulation has been carried out for both CFD and LB algorithms. The Landau-Lifshitz Navier Stokes (LLNS) equations include thermal fluctuations in the stress tensor but numerical schemes to solve them are not nearly as advanced as are the standard CFD solvers [131, 51, 19]. Fluctuations have been included in LB and do not pose any particular numerical problems [3]. Fluctuations have also been included in incompressible solvers in conjunction with the Immersed Boundary Method [12, 83]. The ability to turn fluctuations on or off is an important advantage of continuum-based methods over particle methods. Coupling with the Polymer Chains Regardless of what continuum method is employed, it is necessary to couple that method to the MD description of the polymer chains. The simplest and most com-

54 monly used coupling scheme is to approximate the beads as points and assume for the solvent-induced force on the polymer beads the Stokes-Langevin form F = −6πRH ηvf +FS , where vf is an estimate of the local fluid velocity and FS is an uncorrelated stochastic force whose magnitude obeys the fluctuation-dissipation theorem [146, 58]. This approximation is similar to that in Brownian dynamics, namely, Stokes law is only valid in quasi-static continuum situations, relying on the separations of time and length scales which are usually only marginally separated in realistic situations. Typically the strength of the coupling, RH , is empirically tuned to reproduce experimental measurements. The coupling can also be dealt with when the beads occupy an actual volume, free of fluid. Then stick or slip boundary condition at the surface of the beads are employed, as in both NS [131] and LB [150] simulations of colloidal dispersions. However, these methods are rarely used in polymer simulations due to the complexity when many moving particles are involved, because, the grid size needs to be smaller than the bead size and may need to be adaptively changed when the bead moves. A different alternative is provided by the Immersed Boundary method [12], where the fluid occupies the whole space and the particles, represented as immersed structures, move together with the fluid with a velocity that is a localized average of the fluid velocity. This eliminates the bead inertia from the problem and the need to explicitly enforce boundary conditions on the surface of the beads. The method can be seen as an alternative to Brownian dynamics that correctly captures time-dependent momentum transport in the fluid by explicitly representing the fluid flow, and also includes thermodynamically-consistent thermal fluctuations [83].

55

5.2.3

Explicit Solvent: Particle Methods

An alternative to continuum methods is to use a particle representation of the fluid. The most detailed description is a MD simulation of both the fluid and the solvent. Unlike the classical NS equations, MD automatically and correctly includes fluctuations, internal fluid structure, diffusion, and non-linear transport. Particle methods are also typically simple to implement and can easily accommodate complex boundary conditions.Typically a truncated repulsive Lenard-Jones potential is used for the solvent-solvent interactions. However, even with massive parallelization such MD simulations are limited to short total times and therefore efforts have been made to coarse-grain the solvent to a mesoscopic representation. There, the fluid particles are no longer representative of solvent molecules, but are larger having different dynamics and interactions with each other. However, the viscosity and the stress fluctuations in the solvent must be reproduced correctly. There are mesoscopic particle solvents of progressively decreasing level of microscopic fidelity, and thus increasing efficiency. The handling of the coupling between the solvent and the beads is a separate issue, like for continuum solvents. A particle solvent may be coupled to a polymer chain by including explicit short-ranged solvent-bead continuous [89] or hard-spheres [42] interaction potentials. Efficiency can further be gained by coarse graining the beadsolvent interactions as well, typically using the same ideas as used to coarse grain the solvent-solvent interactions [82, 98]. Dissipative Particle Dynamics (DPD) [111] further coarsens the solvent molecules to obtain a system of weakly-repulsive spheres interacting with a mixture of conservative, stochastic, and dissipative forces. The conservative forces can be used to reproduce the solvent equation of state, while the dissipative forces model viscous friction. The stochastic forces act as a thermostat

56 that ensures detailed balance and correct thermal fluctuations in the DPD fluid. The method has great flexibility and requires significantly less solvent particles and larger time-steps than classical MD, however, it still requires costly integration of differential equations of motion for each of the solvent particles. Such integration of ODEs can be avoided by using a kinetic Monte Carlo method, such as Direct Simulation Monte Carlo (DSMC), to represent the solvent-solvent interactions. The idea is to use stochastic conservative collisions between nearby solvent particles to represent the exchange of momentum and energy. Both multi-particle collisions [82, 98] and binary collisions [42] have been used, as described in Section 5.3.3. The computational efficiency comes at the cost of neglecting the structure of the solvent, as in continuum methods. Recently a new Stochastic Hard-Sphere Dynamics method has been proposed that also uses uncorrelated stochastic binary collisions but still produces a non-trivial fluid structure and a thermodynamically-consistent non-ideal equation of state, similar to those of a DPD fluid [43].

5.2.4

Coupled Methods

Methods that combine several of the techniques described above into a single concurrently coupled simulation can take advantage of their region of validity. Such a simulation may involve several levels each with a different level of microscopic detail. For example, molecular dynamics with complete atomistic detail and realistic potentials may be used for the polymer chain(s) and nearby solvent. The solvent can then be coarse grained to a mesoscopic particle fluid sufficiently far from any chains. The particle method can then be coupled to an explicit fluctuating hydrodynamic description with a fine grid, for example, LB or a fluctuating NS solver. Finally, the

57 hydro grid can be adaptively coarsened in regions even farther from the chain, and a non-fluctuating continuum solver used. This last macroscopic level can use a different method from the fluctuating hydrodynamics level, for example, it could be an incompressible NS solver. Much remains to be done to enable a truly multiscale simulation capable of bridging from microscopic to macroscopic length and time-scales [155, 39].

5.3

Simulation Methods

In this Section we describe in further technical detail the three different techniques we apply to the tethered polymer problem. The majority of the methodology has been previously published so here we only summarize the essential points and cite the relevant works.

5.3.1

Brownian Dynamics

Details of the DNA model and Brownian dynamics simulation method that we use can be found in Refs. [76, 74]. We discretize a double-stranded DNA molecule into a bead-spring chain composed of Nb beads of radius Rb = 77nm (the unit of length, 1 l.u. = 77nm) connected by Ns = Nb − 1 entropic springs. Each spring represents a DNA segment of 4850 base pairs, so that Nb = 11 corresponds to a stained λ-DNA, which has a contour length of 21 µm. In Brownian dynamics, a force balance on this chain leads to a stochastic differential equation for the dynamics of the chain [105],

∆R = [U +

√ D·F ∂ + · D]∆t + 2B · ∆W kB T ∂R

(5.1)

58 where R is the vector containing bead positions, R = {r1 , ..., rN }, U is the unperturbed velocity field at the bead centers, kB is Boltzmann constant, T is absolute temperature, F is the non-hydrodynamic and non-Brownian forces, and D = B · BT is the diffusion tensor. The components of ∆W are obtained from a real-valued Gaussian distribution with mean zero and variance dt. In a unbounded space, the hydrodynamic interactions (HI) enter the chain dynamics through the diffusion tensor, Dij = kB T [(6πηa)−1 Iδij + Ωij ]

(5.2)

where η is the viscosity of the solvent, a is the bead hydrodynamic radius, I is the unit tensor, δij is the Kronecker delta, and Ω is the HI (Stokeslet or Oseen) tensor. Recent work has provided evidence of hydrodynamic coupling to the wall and experimental validation of the use of point-particle (Stokeslet) hydrodynamic interactions (HI) to describe the motion of Brownian particles near a surface [48]. Therefore, it is essential to have wall corrected HI in the simulation to capture the dynamics of a tethered chain correctly. In a bounded space, like near a solid wall, the HI tensor is modified to, Ωij = (1 − δij )ΩOB (ri − rj ) + ΩW (ri − rj )

(5.3)

where ΩOB is the free-space diffusion tensor, and ΩW is the correction which accounts for the no-slip constraint on the wall. The solution for a Stokeslet above a flat plate given by Blake allows us to calculate ΩW exactly [22]. In a square channel or complex geometries, we need to solve this problem numerically with a finite element method to determine ΩW at a grid of points [78]. Based on this description of near-wall HI, Jendrejack et al. [77] predicted that the DNA molecules migrate away from the wall

59 in shear flow, leading to the formation of depletion layers in the near wall region. This prediction has been verified in recent experiments of dilute DNA solutions undergoing pressure-driven flow in microchannels [4, 30]. In different works, Delgado-Buscalioni used a hybrid particle-continuum model method to describe HI [37] and Schroeder et al. used unbounded space HI [127] to study the motion of a tethered chain. We further assume that the chain is ideal (no self-excluded volume interactions between different beads). The entropic springs connecting the beads obey a worm-like chain law Fsij =

kB T |rj − ri | −2 4|rj − ri | rj − ri [(1 − ) −1+ ] , 2bk Nk,s bk Nk,s bk |rj − ri |

(5.4)

where bk is the Kuhn length for DNA and Nk,s is the number of Kuhn lengths per spring. The physical confinement is taken into account through an empirical beadwall repulsive potential of the form −2 3 Uiwall = Awall b−1 k δwall (hi − δwall ) ,

(5.5)

when hi < δwall , where hi represents the perpendicular distance of bead i from the wall, δwall is the cut-off distance. In this work, we choose Awall = 25kB T and δwall = 1/2

bk Nk,s /2 = 0.24 µm. All of the parameters {a, bk , ν} are the same as used in previous work, where it has been shown to successfully reproduce the static and dynamic properties of DNA with contour length 10µm − 126µm [76, 78]. For each parameter set, the sample size is 30 chains unless otherwise specified. All results are presented for DNA at room temperature in a solvent with a viscosity of 1 cP . To study the dynamics of a tethered chain, beads are labeled from 1 to Nb + 1, starting from the tethered point, as illustrated in Fig. 1. The fluid velocity in the

60 flow direction z is a linear function of distance from the wall in the gradient direction x, vz = γx, ˙ where γ˙ is the shear rate, and vx = 0 and vy = 0. Following common experimental practice, the longest relaxation time is calculated by allowing a chain that is initially stretched using a large shear rate to relax to equilibrium. Near equilibrium, the relaxation time is determined by an exponential decay fit the chain extension along the stretch direction, ¯ 2 i = (X ¯ 2 (0) − hX ¯ 2 ieq ) exp(− t ) + hX ¯ 2 ieq . hX τ

(5.6)

An exponential fit to the autocorrelation of the chain extension (relative to equilibrium) parallel to the wall gives similar results. The relaxation time for our λ-DNA is estimated to be 0.59s at room temperature, which is in good agreement with the experimental result of 0.51s [46] after extrapolating the viscosity to 1 cP .

5.3.2

Lattice-Boltzmann

In addition to Brownian Dynamics, we examine the short time correlations of a tethered polymer in a uniform shear flow using a hybrid Lattice Boltzmann (LB) and Molecular Dynamics (MD) code based on the method by Ahlrichs and Dunweg [5]. The Lattice Boltzmann method is a mesoscopic approach to fluid flow calculation and is based on a discrete version of the Boltzmann equation with enough detail to recover hydrodynamic behavior. The LB equation describes the evolution of a singleparticle distribution function, fi (x, t), which is the mass density of particles moving

61 with velocity ei at a time t and position x on a cubic lattice,

fi (x + ei ∆t, t + ∆t) = fi (x, t) +

X j

  Aij fj (x, t) − fjeq (x, t) .

(5.7)

The set of velocities ei is discrete and chosen such that x + ei ∆t always remains a lattice site. The last term describes the collision process in which the distribution function relaxes to a local equilibrium, for which we utilize the BGK (BhatnagarGross-Krook) approximation to the collision operator, Aij = −τ −1 δij , where τ is a relaxation time. The macroscopic hydrodynamic quantities, density ρ, momentum j = ρu, and momentum flux Π, are computed from moments of the particle distribution function, ρ=

X i

fi , j =

X i

fi ei , and Π =

X i

fi ei ⊗ ei .

(5.8)

The equilibrium distribution depends on the macroscopic variables and its form is given by fieq

 u2 ei · u (ei · u)2 − 2 , (x, t) = wi ρ 1 + 2 + cs 2c4s 2cs 

(5.9)

where the weights wi depend on the particle velocity discretization and are determined √ by mass and momentum conservation. The lattice sound speed is cs = ∆x/ 3∆t, where ∆x is the lattice spacing. In this work we solved the distribution function on the standard D3Q19 lattice [150] where the 19 particle velocity components consist of one rest particle, the 6 nearest neighbors in a simple cubic lattice, and the 12 next nearest neighbors in the [110] directions. The corresponding weights are 1/3, 1/18, and 1/36. The LB method avoids the additional mathematical complexities of Navier-Stokes PDE solvers and is straightforward to parallelize efficiently. Using a

62 Chapman-Enskog expansion, the lattice-Boltzmann equation can recover the NavierStokes equations for small Mach and Knudsen numbers, and, within these limits it is second-order accurate in space and time. Compared to the other two methods we apply to the tethered polymer problem, BD and SEDMD, LB is less efficient in this case since it solves for the solvent in the entire domain, even relatively far from the polymer chain. In the LB calculations, the polymer is represented by 25 point particles joined by finitely extendable nonlinear elastic (FENE) springs and interact through a repulsive Lennard-Jones potential among each other and with the walls. Solvent fluctuations are incorporated by adding a stochastic term to the right hand side of the LB equation. This term introduces fluctuations into the momentum flux in a manner that satisfies the fluctuation-dissipation theorem [150]. Coupling between the LB for the solvent and the MD for the solute is achieved through Stokes drag forces and white-noise stochastic forces acting on the monomers. The first monomer in the chain is tethered to the stationary lower wall in a domain having 36, 22, and 24 lattice sites in the streamwise, spanwise, and wall normal directions. The streamwise and spanwise directions are periodic and the bounding upper wall moves with constant velocity, providing the uniform shear.

5.3.3

Stochastic Event-Driven Molecular Dynamics

In addition to Brownian dynamics and Lattice-Boltzmann, we have also applied a purely particle-based method to the tethered polymer problem. The Stochastic EventDriven Molecular Dynamics (SEDMD) algorithm introduced in Ref. [42] combines Event-Driven Molecular Dynamics (EDMD) for the polymer particles with Direct

63 Simulation Monte Carlo (DSMC) [7] for the solvent particles. In SEDMD, the polymers are represented as chains of hard spheres tethered by square wells. The solvent particles are realistically smaller than the beads and are considered as hard spheres that interact with the polymer beads with the usual hard-core repulsion. The algorithm processes true (deterministic, exact) binary collisions between the solvent particles and the beads, without any approximate coupling or stochastic forcing. However, the solvent particles themselves do not directly interact with each other, that is, they can freely pass through each other as for an ideal gas. Deterministic collisions between the solvent particles are replaced with momentum- and energyconserving stochastic collisions between nearby solvent particles. This gives realistic hydrodynamic behavior and fluctuations in the solvent, with tunable viscosity and thermal conductivity, but without internal fluid structure. A recent modification of the DSMC algorithm can be used to achieve a non-ideal equation of state for the stochastic solvent that is thermodynamically-consistent with the density fluctuations [43]. Hard-sphere models of polymer chains have been used in EDMD simulations for some time [139, 107, 101]. These models typically involve, in addition to the usual hard-core exclusion, additional square well interactions to model chain connectivity. Recent studies have used square well attraction to model the effect of solvent quality [104]. Even more complex square well models have been developed for polymers with chemical structure and it has been demonstrated that such models, despite their apparent simplicity, can successfully reproduce the complex packing structures found in polymer aggregation [101, 107]. Here we use the simplest model of a polymer chain, namely, a linear chain of Nb particles tethered by unbreakable bonds. This is

64 similar to the commonly-used freely jointed bead-spring FENE model model used in time-driven MD. The length of the tethers has been chosen to be 1.1Db , where Db is the diameter of the beads. Several particle methods for hydrodynamics have been described in the literature, such as MD [13], dissipative particle dynamics (DPD) [111], and multi-particle collision dynamics (MPCD) [119, 89]. Molecular dynamics is the most accurate model of the fluid structure and dynamics, however, it is very computationally demanding due to the need to integrate equations of motion with small time steps ∆t and calculate interparticle forces at every time step. The key idea behind DSMC is to replace deterministic interactions between the particles with stochastic momentum exchange (collisions) between nearby particles. The standard DSMC [7] algorithm starts with ′

a time step where particles are propagated advectively, ri = ri + vi ∆t, and sorted into a grid of cells. Then, a certain number Ncoll ∼ Γc Nc (Nc − 1)∆t of stochastic conservative collisions are executed between pairs of particles randomly chosen from the Nc particles inside the cell. For mean free paths comparable to the cell size, the grid of cells should be shifted randomly before each collision step to ensure Galilean invariance. The collision rate Γc and the pairwise probability distributions are chosen based on kinetic theory. In SEDMD the polymer chains and the bead-solvent interactions are handled using hard-sphere event-driven molecular dynamics (EDMD) [6, 44, 139, 104] instead of the time-driven MD (TDMD) widely used for continuous potentials. The essential difference between EDMD and TDMD is that EDMD is asynchronous and there is no time step, instead, collisions between hard particles are explicitly predicted and processed at their exact (to numerical precision) time of occurrence. Since particles move

65 along simple trajectories (straight lines) between collisions, the algorithm does not waste any time simulating motion in between events (collisions). SEDMD combines time-driven DSMC with EDMD by splitting the particles between ED particles and TD particles. Roughly speaking, only the polymer beads and the DSMC particles surrounding them are treated asynchronously as in EDMD. The rest of the DSMC particles that are not even inserted into the event queue. Instead, they are handled using a time-driven (TD) algorithm very similar to that used in traditional DSMC. In three dimensions, a very large number of solvent particles is required to fill the simulation domain. The majority of these particles are far from the polymer chain and they are unlikely to significantly impact or be impacted by the motion of the polymer chain. We therefore approximate the behavior of the solvent particles sufficiently far away from any polymer beads with that of a quasi-equilibrium ensemble. In this ensemble the positions of the particles are as in equilibrium and the velocities follow a local Maxwellian distribution whose mean is the macroscopic local velocity. These particles are not simulated explicitly, rather, we can think of the polymer chain and the surrounding DSMC fluid as being embedded into an infinite reservoir of DSMC particles which enter and leave the simulation domain following the appropriate distributions. Using such open or stochastic boundary conditions dramatically improves the speed, at the cost of small errors due to truncation of hydrodynamic fields. This truncation can be avoided by coupling DSMC to a continuum fluctuating hydrodynamic solver [155]. We have made several runs for different polymer lengths and also bead sizes. One set of runs used either Nb = 25 or 50 large beads each about 10 times larger than a solvent particle. Another set of runs used either Nb = 30 or 60 small beads each

66 identical to a solvent particle, with faster execution but nearly identical results. In the simulations reported here we have used rough wall BCs for collisions between DSMC and non-DSMC particles [42]. This emulates a non-stick boundary condition at the surface of the polymer beads. Using specular (slip) conditions lowers the friction coefficient, but does not qualitatively affect the behavior of tethered polymers. All of the runs used open boundary conditions, where about 153 DSMC cells around each bead were explicitly simulated. Note that for (partially) collapsed polymer chains the total number of explicitly simulated cells is much smaller than 153 Nb . The Nb = 30 runs were run for about 6000τ0 relaxation times, and such a run takes about 6 days on a single 2.4GHz Dual-Core AMD Opteron processor. Even for such long runs the statistical errors due to the strong fluctuations in the polymer conformations are large, especially for correlation functions at long time lags t > τ .

5.4

Results

The main goal of our paper is to reinvestigate the tethered chain problem through extensive long time simulations (thousands of longest relaxation time of the tethered polymer, τ ) involving different representations of polymer and solvent, including Brownian dynamics [74] (BD), the Lattice Boltzmann method [5] (LBM), and a recent Stochastic Event-Driven Molecular Dynamics [42] (SEDMD) algorithm. In this section we present comparison results from our simulations. More extensive results for the tethered polymer problem obtained using the SEDMD algorithm are presented in Ref. [42]. Since the three different methods that we use give similar results and Brownian Dynamics is the fastest methodology, the majority of the results we present

67 will be from BD simulations with Nb = 11 (Ns = 10), unless otherwise indicated. Of the three methods used here, LB is the slowest and thus the LB results are of more limited duration. We emphasize that direct computational comparison between the methods is unfair. Most significantly, the LB runs use periodic boundary conditions and have to fill the whole simulation domain with explicit solvent (lattice points). By contrast, the SEDMD runs use open boundaries and thus use much less explicit solvent, whereas the Brownian dynamics does not use an explicit solvent at all. Doyle et al. proposed a cyclic dynamics mechanism for a tethered polymer chain in shear flow (Fig. 5.1) based on Brownian dynamics simulation results [46]. According to this scenario, when thermal fluctuations cause motion in the gradient direction x (from state 1 to state 2), the chain is driven away from the wall and experiences higher hydrodynamic drag. This leads to further stretching and an increase of the extension in the flow direction z (state 3). Due to the finite extensibility of the chain, the extension in the z direction is finite and depends on the shear rate and chain properties. After stretching, the coupled torque of the hydrodynamic drag and spring forces will rotate the chain towards the wall (state 4). As the chain get closer to the wall, the flow velocity decreases and entropic recoiling becomes dominant, resulting in a decrease of the z extension (state 1). The tethered chain could take other dynamical paths than following the one described above, such as restretching or recoiling after state 2 by random motion in −x and −z direction, respectively. In Fig. 5.2 we show the probability distribution function (pdf) ρ(z, x) of the end bead in the z − x plane at Wi = 0 and Wi = 2 for the three different methods. The results are presented in dimensionless units by normalizing the unit of length by the average radius of gyration in the x direction at Wi = 0. Here, again, we want to

68

4.0

4.0

4.0

3.5

3.5

3.5

3.0

3.0

3.0

2.5

2.5

2.5

2.0

2.0

2.0

1.5

1.5

1.5

1.0

1.0

1.0

0.5

0.5

0.5

0.0

0.0 -3

-2

-1

0

1

2

3

0 0.033 0.065 0.098 0.13 0.16 0.20 0.23 0.26

0.0 -3

-2

-1

0

1

2

3

-3

-2

-1

0

1

2

3

(a) Wi = 0 3.0

3.0

3.0

2.5

2.5

2.5

2.0

2.0

2.0

1.5

1.5

1.5

1.0

1.0

1.0

0.5

0.5

0.5

0.0

0.0 -1

0

1

2

3

4

0 0.081 0.16 0.24 0.33 0.41 0.49 0.57 0.65

0.0 -1

0

1

2

3

4

-1

0

1

2

3

4

( b ) Wi = 2 Figure 5.2: Probability distribution of the end bead of the tethered DNA molecule in a dimensionless x − z plane at Wi = 0 and Wi = 2. The visible differences can likely be attributed to the differences in the boundary conditions between the different methods, as well as the different elasticity of the chains. (Left) Brownian dynamics. (Middle) Stochastic Event-Driven Molecular Dynamics. (Right) Lattice Boltzmann Method.

69 emphasize that we are not expecting perfect match between methods. In particular, the different methods implement different effective boundary conditions at the wall surface. In Brownian dynamics, an essentially reflective boundary condition appears, while in the case of a hard-sphere chain a perfectly reflective boundary condition is appropriate. For the LB runs an intermediate case appears, where the repulsion from the wall is stronger than hard spheres but still finite-ranged. These boundary effects are clearly visible in the results in Fig. 5.2, where the BD results show a depletion layer near the wall where as the SEDMD and LB results show the bead spending more time near the wall. In Fig. 5.3 we compare the dependence of the relaxation times τx/y/z along the three different axes on the flow rate among the three different methods. The figure shows reasonable agreement between the different techniques, especially considering the large errors inherent in determining relaxation times. We calculate the relaxation times by fitting an exponential decay to the intermediate portion of the autocorrelation function 0.2 < C(t) < 0.8 of the position of the end bead along the three coordinate axes. The LB calculations use periodic boundaries with a narrower box in the spanwise (y) direction than in the streamwise (x) direction, which makes the relaxation times τx (Wi = 0) and τy (Wi = 0) unequal, as they must be by symmetry. We have scaled τy (Wi) (the shorter axes) in the LB results by a constant factor so as to correct this strong boundary effect at Wi = 0. Among the three relaxation times, the relaxation in the direction perpendicular to the wall τz is the shortest, even for no flow. Note than in this work, following Ref. [42], the relaxation time along the flow direction τx is used to define the internal relaxation time and thus Wi when comparing among the different methods. Note also that it is τx that seems to get

70

circles=BD, squares=SEDMD, diamonds=LB 2 2 x z y

1.5

τ z / τx

1.5

τ y / τx

τ / τ0

1 0.5 1 0

0

1

2

3

4

5

6

0.5

0

0

1

2

3

4

5

6

Wi

Figure 5.3: Dependence of the dimensionless relaxation time τ (Wi)/τ (Wi = 0) of the tethered chain along the three coordinate axes as a function of dimensionless flow rate Wi. The inset shows the ratios of the different relaxation times as a function of Wi.

71 most strongly reduced as Wi increases. To study the time scale associated with the fluctuating process (cycle) quantitatively and to find the correlation between different chain segments, we calculated the cross-correlation functions (CCF) of beads’ positions. We also calculated the power spectral density (PSD) in search of periodicity. The CCF and PSD are the natural tools for examining the relationship between two time dependent random variables in the time and frequency domain respectively. The mean-removed CCF of two time series α(t) and β(t) is defined as

Cαβ (T ) =

¯ E[(α(t + T ) − α)(β(t) ¯ − β)] σα σβ

(5.10)

where α ¯ = E(α) is the mean, σα2 = E(α2 )−[E(α)]2 is the standard deviation, and T is the time lag. A significant peak in the CCF at lag T indicates that α(t) is correlated to β(t) when delayed by time T . In the frequency domain, the PSD is the norm of the Fourier transform of the CCF,

Z

Sαβ (ν) =

∞

−∞

Cαβ (T )exp(−2iπνT )dT

(5.11)

Note that this is the standard definition used in the engineering literature, and here the frequency ν = 1/T is actual frequency (inverse period) rather than angular frequency ω = 2πν. To produce a PSD with accurate sampling around interesting frequencies, long simulation times and a high sampling frequency are essential. We examined various choices of variables to represent the motion of the chain and have found little qualitative difference between them. We have chosen the position of the end bead rNb = (x, y, z) as be the best option [42]. Extensive computational efforts

72 have been undertaken to determine the CCF and PSD of the end bead coordinates as function of chain length N and shear flow parameter Wi. The CCF Czx (t) of the end bead at various Wi is shown in Fig. 5.4(a). The shape of the CCF is consistent with the cyclic dynamics mechanism proposed by Doyle. Clearly, in the absence of flow, Wi = 0, the movements in the x and z directions are uncorrelated on all time scales. When shear flow is introduced, the movements in flow direction and gradient direction are coupled together due to the nature of the flow and the finite extensibility of the chain, as reflected in the rise of a prominent peak in the CCF. When thermal motion in +x direction occurs, the chain will be stretched with an increase in +z, which leads to a positive correlation. Similarly, when motion in −x direction is introduced, the chain will recoil in the −z direction as the drag decreases, which also leads to a positive correlation. As expected, the larger the shear rate, the greater the correlation. There’s only one significant peak in the long time correlation function, shown in the inlet of Fig. 5.4(a), which suggests that all correlations are short-lived and not periodic. Turning attention to the correlation between different chain segments, Fig. 5.5 shows the CCFs of several beads along the chain at Wi = 2. One striking feature is that for beads sufficiently far from the tether all curves pass the time axis at the same time lag. The fact that all CCFs have the same shape indicates that a common movement pattern exists for the whole chain. The inset in Fig. 5.5(a) shows the CCFs of the x coordinates of different beads. Although the correlation decays as the distance along the chain increases, it confirms that all beads move in a cooperative manner, indicated by the fact that the peak positions are all at zero time lag. The CCF for the end bead for various chain lengths are compared in Fig. 5.5(b), to show that there is no fundamental difference

73

Figure 5.4: (a) Normalized cross-correlation functions (CCF) Czx (t) of the end bead’s coordinates in flow direction z and gradient direction x as a function of nondimensional time, at various Wi for Ns = 10. The inset shows longer time lags. (b) Power spectral density (PSD) Szx (ν) of the end bead’s coordinates as a function of non-dimensional frequency. The results are averaged over 30 runs for a total simulation time is 103 τ , and 104 τ for Wi = 5.

74

Figure 5.5: (a). CCFs of end beads’ coordinates at Wi = 2 for a chain with Ns = 10 and simulation time is 1000 τ . The number in the legend is the bead label as shown in Fig. 5.1. The inset shows the CCFs for x coordinates of different beads to study the correlation of the dynamics between different beads (similar results are obtained for the y axes). (b) CCFs of end bead as function of chain length at Wi = 5. The number in the legend is the number of springs Ns in the chain. The inset shows longer time lags.

75 between different chain lengths, ranging from 20 µm (Ns = 10) to 80 µm (Ns = 40). We have also established that these results are insensitive to the cut-off distance and the magnitude of the repulsive potential between the wall and chain segments. In Fig. 5.6 we compare the cross-correlation function Czx (t) at Wi = 2 among the three different methods: Brownian Dynamics, Lattice-Boltzmann, and Stochastic Event-Driven Molecular Dynamics. In particular, our goal is to verify the pervasive assumption that the dynamics of polymer chains in shear flow is essentially universally quantitatively determined by Wi for a wide range of flexible chains. Furthermore, it is important to cross-validate the different methods against each other, given that each of them makes certain assumptions and has somewhat different range of applicability. The results in Fig. 5.6 indeed show reasonable agreement between the different methods. Perfect agreement is not expected because the polymer models are different among the different methods. The cross-correlations we measure are not consistent with periodic motion. The PSD calculation does not show discernible peaks either, as shown in Fig. 5.4(b). All that we can reliably extract from the results is that the response of the chain to a large thermal fluctuation (the “cycle”) is reproducible for short times, and we find no evidence of sustained correlations (oscillations) at times longer than the internal relaxation time of the chain. For a free chain in shear flow, where rotations of the chain are possible, one can count the number of tumbling events per unit time and define that as a cycling time. The distribution of the delays between successive tumbling events is itself important. If this distribution is sharply peaked, that would be consistent with a periodic motion with a well-defined period. If the distribution is exponential, this would indicate a Poisson-like tumbling process. Several recent works have proposed an exponential distribution for the delay between

76

C zx

0.1

BD Wi=2.0 (N=20, error~0.005) SEDMD Wi=1.8 (N=30, error~0.01) LB Wi=2.0 (N=25, error~0.05)

0

-0.1 -2.5

-2

-1.5

-1

0

-0.5

0.5

1

1.5

2

t/τ

Figure 5.6: Comparison of the cross-correlation function Czx (t) at Wi = 2 among the three different methods: Brownian Dynamics (30 runs about ∼ 1000τ long), LatticeBoltzmann (run is ∼ 600τ long), and Stochastic Event-Driven Molecular Dynamics (10 runs about ∼ 1000τ long).

77 successive tumblings [57, 34, 156]. Furthermore, the tumbling time was found to be related to the internal relaxation time of the chain [34, 156]. For tethered chains, we cannot even identify and count a unique event such as tumbling and thus we cannot extract a repetition frequency for the “cycle”. In Appendix C, we analytically calculate the CCF for a Brownian particle tethered to the origin with a harmonic spring and subjected to shear flow. This simple dimer model qualitatively reproduces the features we see in the CCF for the tethered chains, namely, a single peak at t ∼ τ of width ∼ τ and height ∼ Wi. Better quantitative agreement is obtained when a nonlinear spring and a hard wall surface are also included (without hydrodynamics). The PSD for the dimer model shows no peaks and there is only a single time-scale in the dynamics, namely, the intrinsic relaxation time τ . Furthermore, the analytical form of the CCF shows that by a slight modification of a tunable parameter one can obtain a CCF fully-consistent with our numerical results for longer chains. This analytical CCF has an analytical PSD that does show a broad peak at small frequencies ντ ∼ 0.1, very similar to the previously reported peaks used to justify the claims to periodicity in the chain motion [127, 37, 38]. This peak is weak and broad even when plotted on a logarithmic axes and its exact shape and maximum will vary depending on the particular model, variables used in calculating the PSD, Wi, the definition used for calculating τ and Wi, etc. We therefore believe that its interpretation as evidence of periodic motion is not justified. The calculations in Appendix C for a dimer in shear flow also demonstrate that a qualitatively similar behavior is observed even without hydrodynamic interactions. Our results from Brownian Dynamics simulations in the free-draining limit confirm this and show that the HI do not affect the results significantly, so long as the relax-

78 ation time is recalculated when computing Wi. In the tethered case, we believe that the competition between frictional and elastic restoring forcing dominates and the hydrodynamic interactions are a weak perturbation. Therefore, it is not surprising that the proper inclusion of hydrodynamic interactions is not essential for the tethered polymer problem, as reasoned theoretically for a free chain in shear flow in Ref. [156].

5.5

Conclusions

We studied the dynamics of a polymer molecule tethered to a hard wall and subjected to a shear flow. We found consistent results among three methods utilizing different representations of the solvent, Brownian Dynamics (BD), Lattice-Boltzmann (LB), and Stochastic Event-Driven Molecular Dynamics (SEDMD). Specifically, BD implicitly represents the solvent, LB explicitly represents the solvent flow on a discrete lattice, and SEDMD utilizes a particle-based solvent. The three methods also utilized different polymer chains, namely, the BD simulations used a worm-like chain, the LB simulations used a FENE-LJ chain, and for SEDMD we used a tethered chain of hard spheres. The correlation functions of the position of the end bead question the existence of periodic motion, as previously suggested. The cross-correlation function between the bead positions along the flow and gradient directions shows a single peak indicative of a fluctuation-dissipation cycle of duration comparable to the relaxation time of the polymer. The corresponding Fourier representation, the power-spectral density, shows no peaks. We find that neither the chain length of the polymer N, nor the

79 dimensionless shear rate Wi, qualitatively alter the results, and in the Appendix C we give some calculations for a very simple model of a dimer in shear flow that reproduces the essential features of the observed peak in the cross-correlation function. While our conclusions are rather different from other authors, our results are statistically consistent with those presented in the literature. Specifically, the shape and position of the peaks in the cross-correlation functions are very similar to reported results, however, we did not observe large oscillations in the CCFs previously identified as signatures of periodic motion [37]. We believe that this is due to the requirement of very long simulation times to obtain good statistics for the time-correlation functions at long time lags, as necessary to establish periodicity. Not all previous studies have been able to reach sufficiently long simulation times. Another important point we clarified is that maxima in the power-spectral density does not necessarily indicates a periodic motion, which we demonstrate in Section C using an analytic dimer model. Namely, an analytical shape is suggested by the dimer calculations that can exhibit peaks very similar to those reported in the literature through small adjustments of a tunable parameter, whose appropriate value likely depends on details of the model used and the exact variables used in the calculations of the power spectrum. Furthermore, different if not conflicting ways have been used to define and calculate the “cycling time”, without properly distinguishing between the duration of a cycle and the interval between cycles. Even more importantly, the very concept of a cycle in the chain motion as a well-defined countable event, analogous to the case of a free chain in shear flow, should be questioned. Our results are consistent with a simple traditional picture of continuous thermal fluctuations dissipated by deterministic friction, leading to exponentially-decaying correlation functions.

80

Chapter 6 Flow-driven DNA dynamics in complex geometries In this chapter, we consider the dynamics of a flow-driven DNA molecule in micro/nanofluidic devices. In this case, the configuration-dependence of the mobility tensor cannot be ignored and the solvent velocity field is in general non-linear on the length scale of the molecule. We present an immersed boundary method that allows fast Brownian dynamics simulation of polymer chains and other particles in complex geometries with fluctuating hydrodynamics.1 This approach is applied to study the dynamics of a flow-driven DNA molecule through a nanofluidic slit with an embedded array of nanopits. We investigate the dynamics of the DNA molecule as a function of the Peclet number and chain length, as well as the influence of hydrodynamic interactions by comparing with free draining simulation results. 1

This chapter is based on the manuscript: Y. Zhang, J. J. de Pablo, and M. D. Graham, submitted to J. Chem. Phys. (April 12, 2011)

81

6.1

Introduction

Transport of polymer solutions in constricted spaces [1] is a long-standing research topic with many applications including polymer enhanced-oil-recovery, size exclusion chromatography [153] and gel-electrophoresis [151] and, recently, single DNA molecule analysis using micro- and nano-fluidic devices [140, 125]. In many fluidic devices designed for separation, manipulation, and sequencing of DNA, the critical dimension of the constriction approaches the radius of gyration of the polymer molecule or smaller. In particular, capturing the interaction between polymer dynamics and fluid motion in a confined geometry is essential to a proper description of these devices.[1] The purpose of the present work is to introduce a computational approach to the dynamics of polymers in complex confined geometries, and to apply the approach to an interesting recent set of experiments on the flow of DNA solutions over arrays of nanopits [36]. Many devices with well-defined microstructure have been proposed for DNA manipulation, especially for separation and stretching purposes using electrophoresis[45]. By contrast, much less attention has been paid on pressure (flow) driven DNA dynamics in microfabricated devices [141, 36, 143]. There are significant differences between the electrophoresis case and the pressure driven case. Consider for example DNA through a slit geometry. In electrophoresis, in a uniform field, the velocity of each DNA segment is the same everywhere in the channel; by contrast, the unperturbed velocity profile is parabolic across the channel in the pressure driven case at low Reynolds number, and this velocity gradient can give rise to interesting transport phenomena, such as Taylor dispersion [141]. Furthermore, hydrodynamic interactions (HI) between objects such as polymer segments in an unconfined domain are long-

82 ranged, leading to strong many-body effects, for example, Zimm scaling of the selfdiffusion coefficient of a polymer chain. In confined geometries, the long-ranged nature of hydrodynamic interactions changes substantially, leading to significant changes in polymer dynamics. [1]. In a slit geometry, for example, the velocity field due to a point force perpendicular to the walls decays exponentially. For a force parallel to the walls, the velocity field has a parabolic form in the wall-normal direction and decays as 1/r 2. This is still a slow decay, but the symmetry of the flow leads to cancellations upon averaging that result in screening [9, 144, 14]. Both experiments and simulations of the diffusion of long flexible DNA molecules in slits [31, 78, 75] are consistent with screening of hydrodynamic interactions on the scale of slit height. Nevertheless, on smaller scales, HI are not screened and lead to changes in segment mobilities and boundary effects such as cross-stream migration, which leads to depletion layers much larger than the equilibrium chain size [77, 90, 30, 68, 29, 81]. Several recent computational studies investigated the effect of HI on dynamics of polymer molecules in complex geometries.[54, 73, 93, 66, 154, 29, 56] Among the problems studied in these works, forced translocation of a polymer through a nanopore has received substantial attention because of its potential applications to rapid DNA sequencing. In this problem of translocation, one end of a polymer molecule is placed at the entrance to a nanopore embedded in a membrane, and then the polymer molecule is pulled through the pore by a constant force exerted on that leading end segment or on the chain segments in the pore. Using different simulation methods and/or different polymer models, Izmitli et al. [73], Hernandez-Ortiz et al.[66], and Fyta et al.[54] all found that HI shortens the translocation time and alters the scaling exponents for the power-law dependence of the translocation time on polymer length

83 [54, 73]. Using the general geometry Ewald-like method (GGEM)[67], HernandezOrtiz et al. studied the effect of HI on polymer solutions with finite concentration flowing over a channel with grooves oriented perpendicular to the flow direction was considered.[69] At low concentration and high Weissenberg number W i, which is defined as the polymer relaxation time times fluid deformation rate, the groove was almost completely depleted of polymer chains. At concentration approaching overlap, the concentration difference between the bulk and groove is substantially reduced, but only if hydrodynamic interactions were included in the simulation. All these studies suggest that HI is not only important but essential in the simulation of polymer solutions in complex geometries. The specific system of interest in the present work is a nanoslit device that exploits partial confinement to manipulate a DNA molecule. [116, 36] In the experiments, the height of the channel is H = 100 nm, and the bottom wall of the slit contains square nanopits (100 nm deep and 500–1000 nm wide). When a solution of λ−DNA, which has a radius of gyration Rg of 770 nm, is introduced into such a system, the DNA is highly confined as reflected by the confinement ratio Rg /H = 7.7. Because the confinement is weaker in the regions of the slit where there is a pit, DNA segments are preferentially found there. This entropic well tends to pin sections of the chain in these regions (Fig. 6.1(b)). When multiple pits are present (and not too far apart, 1-2 µm in experiments), stretched sections of DNA will span the more confined regions between pairs of pits. When a pressure drop is applied across the channel, the DNA molecule is transported downstream by the flow, hopping between neighboring pits (Fig. 6.1(c)). When a molecule gets trapped in a pit, a free end of the molecule first diffuses out of the pit, then is dragged downstream by the flow. When the free

84 end hits the next pit, the rest of the chain is dragged along into the second pit. In Ref. [36], an equilibrium model was proposed, based on the hypothesis that the dynamics of DNA at low Peclet number, which is defined as the ratio between the diffusive and convective time scale, is thermally activated, and thus the Kramer’s theory of barrier crossing could be applied. There are two important issues that were not addressed in Ref. [36]. First, no direct evidence was shown to support the hypothesis that the dynamics is thermally activated. Specifically, no analysis about the distribution of the residence time in each pit as a function of Peclet number (fluid velocity) was performed. Second, how the chain dynamics change with DNA size was not systematically examined. We address these issues using an accelerated Brownian dynamics method for modeling the hydrodynamics of polymer solutions in complex geometries. This method is based on Green’s function solutions of the Stokes equations, combined with a fast Stokes solver and a variant of the immersed boundary method. This paper is organized as follows: A general discussion of the wide range of techniques for modeling the hydrodynamics of polymer solutions in complex geometries is given in Sec. 6.2. In Sec. 6.3, we present the DNA model, governing equations, and the new immersed boundary method that we have developed to calculate fluctuating hydrodynamics in complex geometries. Sec. 6.4 gives results from the application of the algorithm to the problem of flow-driven DNA dynamics in a nanofluidic device with embedded nanopits array, and some concluding remarks are given in Sec. 6.5.

85

Figure 6.1: DNA flowing through a nanoslit with embedded nanopit arrays (reproduced with permission from Ref. [36]. Copyright 2009 by the Institute of Physics) (a) Scanning electron micrograph of typical 1 × 1 µm square pits embedded in a nanoslit at 1 µm intervals. (b) Epifluorescent images of stained λ-DNA molecules confined in a 107 nm deep slit embedded with 1 × 1 µm pits spaced 2 µm apart. (c) Fluorescent images of λ-DNA travelling across a nanopit array under an applied pressure of 40 mbar in the same device as shown in (b).

86

6.2

Methods for hydrodynamics of confined polymer solutions

In this section we give a brief overview of various methodologies for modeling hydrodynamics of polymer solutions in complex geometries. For theoretical analysis, the reader is referred to Graham.[1] The various methods for computational hydrodynamics of polymer solutions can be divided into two major categories. The first category is purely continuum methods that use constitutive equations for the polymer solution. These models only apply at macroscopic scales, when the number of polymer chains in an elementary fluid flow volume is large, so that statistical averages of the chain conformations can be used as parameters in constitutive models of the time-dependent stress as a function of the strain rate history. The second major category of methods explicitly simulates the motion of each polymer chain using a coarse-grained (bead-spring) model of the chains. Multibead representations of the chains are capable of complex chain conformations but require models for the bead-bead interactions. There are also two major types of approaches for dealing with the fluid dynamics of the solvent and thus the hydrodynamic interactions between chain segments. One approach, exemplified by Brownian dynamics (BD) and Lattice Boltzmann (LB) methods, uses a continuum model (Stokes, Navier-Stokes or Boltzmann equation) to represent the solvent. The other uses a discrete, particle-level, model to represent the solvent; examples include full atomistic molecular dynamics (MD), Stochastic Rotational Dynamics (SRD), and Dissipative Particle Dynamics (DPD). As we are interested in long-time (¿ 1 second) dynamics of a single polymer molecule in complex geometries, this brief review focuses

87 on mesoscopic methods, namely, BD, LB, SRD, and DPD, especially their different treatments of hydrodynamics, boundary conditions, and Brownian fluctuations. In Brownian dynamics simulations of polymer molecules in unbounded or periodic domain, the HI are usually incorporated using analytical Green’s functions for Stokes’ equations: I.e. the beads exert point forces on the fluid and Green’s functions are used to determine the fluid velocity at any other bead position. For a system with N beads, these interactions are incorporated in the N ×N mobility matrix[105]. While for true point forces, the Green’s function for the Stokes equations – the Oseen tensor – is appropriate, in Brownian dynamics simulation it is necessary to use regularized point forces. The most common regularization is the Rotne-Prager-Yamakawa tensor [120, 157], while in the present work a different, computationally convenient regularization is used. Corresponding to the mobility matrix representation of the HI, a matrix formulation of the fluctuation-dissipation theorem – a factorization of the mobility tensor – must be implemented. Ermak and McCammon suggested a Cholesky factorization which requires O(N 3 ) operations[49], while Fixman [53] observed that the square root matrix of the mobility matrix can be accurately approximated via Chebyshev polynomials, and this algorithm brings down the computation complexity to roughly O(N 2 ) or O(N 2.25 ).[76] Another advantage of Fixman’s algorithm, one that is essential to this study, is that it works even when an analytical form of the Green’s function is not available, as is the case for flow in complex geometries. For periodic domains, the scaling can be dramatically improved, to O(NlogN), by using particle-mesh Ewald (PME) methods, which are based on Hasimoto’s solution for the Green’s function for Stokes flow driven by a periodic array of point forces [65, 16]. BD can also be efficiently applied to nonperiodic geometries. In simple cases

88 like flow in a half-plane bounded by a no-slip wall, the exact Green’s function for that geometry can be used. To handle more complex boundary conditions, such as indentations or bumps on walls, it requires incorporating other methods such as finite element methods [78, 69], which can be done in the framework of the general geometry Ewald-like method (GGEM) [67], which also has O(N log N) or O(N) scaling. In the present work, we combine GGEM with a variant of the Immersed Boundary Method.[110, 97, 114] The lattice Boltzmann method is a mesoscopic approach to fluid flow calculation and is based on a discrete version of the Boltzmann equation with enough detail to recover hydrodynamic behavior [47]. The LBM avoids the some of the computational complexities of Navier-Stokes (NS) solvers and is straightforward to parallelize efficiently. Using a Chapman–Enskog expansion, the LB equation can recover the NS equations for small finite Mach and Knudsen numbers, and, within these limits it is second-order accurate in space and time. Coupling between the LB for the solvent and the molecular dynamics for the polymer is achieved through Stokes drag forces. The Brownian fluctuation is introduced via random fluctuations added to the stress tensor and the beads dynamics directly. [47] Usually, the no-slip boundary is introduced by a bounce-back collision rule in which incoming fluid particles are reflected back towards the nodes from which they originated, and this results in no-slip boundary conditions. For complex boundaries, however, the boundary implementation may not be trivial [85, 33, 28]. The Lattice Boltzmann-molecular dynamics scheme has been successfully applied to the problems of polymer translocation through a pore [54, 73, 93], polymer cyclic dynamics near a solid wall [158], and polymer migration in a square channel [81].

89 Stochastic rotation dynamics is a relatively recent particle-based mesoscopic simulation method. In the SRD algorithm, the fluid is modelled as a system of point particles. The algorithm consists of two steps, the streaming step and the collision step.[80, 59] In the streaming step, the particles move ballistically during a certain amount of time. In the collision step, the particles are sorted into cubic cells of a certain size. In each cell, the velocity vector of every particle is rotated an angle α about a randomly chosen axis, relative to the center-of-mass velocity of all the particles within that particular cell. In addition, a random shift is added in order to ensure Galilean invariance. The polymer segments are coupled to the fluid by including them in the collision step. The streaming and collision steps are designed to conserve mass, momentum, and energy. In the long time limit, correct hydrodynamics are recovered. SRD also naturally contains Brownian fluctuations. To implement no-slip boundary conditions in SRD, there are two approaches: the bounce-back scheme and the stochastic scheme. The bounce-back scheme was first proposed by Malevanets and Kapral [91], and is a direct analog of the no-slip boundary condition commonly applied in LBM. Lamura and Gompper found that the simple bounce-back rule is not sufficient to guarantee no-slip without corrections. [87] In the stochastic scheme, proposed by Inoue et al., the velocity of the SRD particle that comes into a wall cell is changed to a random velocity obeying a Maxwell-Boltzmann distribution[72]. This boundary condition, however, allows flow to penetrate small objects.[60] Watari et al.[154] and Chelakkot et al.[29] have used this method to study the dynamics of polymer solutions between two parallel planes, and qualitatively reproduced the hydrodynamic migration phenomenon. Dissipative particle dynamics coarse-grains the solvent molecules to obtain a sys-

90 tem of weakly repulsive (soft) spheres interacting with a mixture of conservative, stochastic, and dissipative forces.[50, 64, 118] The conservative forces can be used to reproduce the solvent equation of state, while the dissipative forces model viscous friction. The stochastic forces act as a thermostat that ensures detailed balance and correct thermal fluctuations in the DPD fluid. The method has great flexibility and requires significantly fewer solvent particles and larger time steps than classical MD. The influence of solid boundaries is typically very strong in DPD simulations. The modeling of no-slip boundaries in the context of DPD has been addressed in the past by freezing the particles that represent the solid boundaries[117], and by modifying the repulsive forces[79]. While these efforts introduce minimal additional algorithmic complexity to the standard DPD formulation, care has to be taken to avoid numerical artifacts such as depletion of particles near the wall or artificial ordering of the near-wall particles leading to large density fluctuations[112]. To study the dynamics of polymer solutions in complex geometries, the ability to treat no-slip boundary conditions efficiently and accurately becomes important. In all three mesoscopic models (LBM, SRD, and DPD), it is not a trivial task to correctly implement the no-slip boundary conditions. Another difference between BD and the three mesoscopic models is that those methods describe compressible inertial fluids in which hydrodynamics need time to propagate through the solution; the Reynolds and Mach numbers Re and Ma are finite. In Green’s function-based methods, on the other hand, the condition Re ≪ 1, Ma ≪ 1 are satisfied explicitly through the use of the Stokes equation. In the present study, we propose an immersed boundary method-based Green’s function method, where the boundary is represented by many regularized points and

91 the Stokes equation is solved on a Cartesian grid. The validity and the accuracy of the new method are scrutinized by simulating laminar flow in a slit, uniform flow around a solid sphere, and flow generated by a Stokeslet near a single wall.

6.3 6.3.1

Polymer model and simulation method Model and governing equations

In this section, we present the coarse-grained model used in this work to study the dynamics of polymer in complex geometries. We focus attention on a model of λDNA[158, 76, 74], described by a bead-spring chain model composed of Nb beads of hydrodynamic radius a = 77 nm connected by Ns = Nb − 1 entropic springs. a is also the length unit used through out this work for nondimensionlization. Each bead represents a DNA segment of 4850 base pairs, i.e., Nb = 11 corresponds to a fluorescently stained λ-DNA, which has a contour length of 21 µm and radius of gyration of 770 nm.[136, 137] The length unit in this work is a = 77 nm. The spring connecting the beads i, j obeys a worm-like chain force law [94] fijs =

|rj − ri | −2 4|rj − ri | rj − ri kB T [(1 − ) −1+ ] . 2bk Nk,s bk Nk,s bk |rj − ri |

(6.1)

Here, bk is the Kuhn length for DNA and Nk,s is the number of Kuhn length per spring. The physical confinement, or the steric interaction between DNA and wall, is

92 taken into account through an empirical bead-wall repulsive potential of the form

Uiwall

=

   A

−1 −2 wall bk δwall (hi

− δwall )3 , hi < δwall

  0,

(6.2)

hi > δwall.

Here hi represents the perpendicular distance of bead i from the wall, δwall is the cutoff distance. In this work, we choose Awall = 50kB T /3 and δwall = 1a. The excluded volume potential between two distinct beads is −3|ri − rj |2 3 3/2 1 2 ) exp( ), ( Uijev = νkB T Nk,s 2 4πSs2 4Ss2

(6.3)

where ν is the excluded volume parameter, and Ss2 = Nk,s b2k /6. All of the parameters {a, bk , ν} are the same as used in previous works[76, 74], where it has been shown to successfully reproduce the static (radius of gyration, Rg ) and dynamic properties (diffusivity D and longest relaxation time τ0 ) of DNA with contour length of 10 − 126 µm.

6.3.2

Governing equations

We are concerned with the numerical simulation of the Brownian motion of a single DNA molecule immersed in an isothermal incompressible Newtonian solvent in the limit of zero Reynolds number. When the DNA molecule is represented by the beadspring model as described above, the equation of motion for each bead is determined by the balance of Stokes drag on the bead due to the surrounding fluid, the Brownian forces exerted by the fluid on the bead, and other non-hydrodynamic, non-Brownian intra/inter molecular forces. The force balance on the beads leads to a stochastic

93 differential equation which governs the dynamics of the chain and can be solved numerically[105]: dR = (U∞ + M · F + kB T

√ ∂ · M)dt + 2B · dW, ∂R

B · BT = kB T M.

(6.4) (6.5)

Here R is the vector containing all bead positions {ri }, kB is the Boltzmann constant and T is the absolute temperature. The vector containing total non-Brownian, nonhydrodynamic forces acting on the beads is denoted F. From Eq. 6.4 we can see that there are four contributions to the change of chain configuration after a time step dt: unperturbed flow velocity U∞ in complex geometries, perturbed velocity M · F due to the polymer forces, drift due to the configuration dependence of the mobility √ ∂ · M, and Brownian displacement 2B · dW. The unperturbed flow velocity tensor ∂R U∞ at R is the velocity when DNA is absent in the system. As noted, the second term in Eq. 6.4, M · F, is the velocity generated by the motion of the polymer beads. It arises from the fact that beads immersed in a fluid generate flows as they move due to various forces, and similarly they move in response to fluid motion through the Stokes drag. As described in last section, we treat bead i of the bead-spring chain as a sphere of hydrodynamic radius a. Then the relationship between the bead velocity ui and the drag force it exerts on the fluid is given by the Stokes law fi = ζ(ui − u(ri )), where ζ is the bead friction coefficient ζ = 6πηa, u(ri ) is the fluid velocity at the bead position, and η is the fluid viscosity. (Note that the finite size of the bead only arises in the friction coefficient.) Through these hydrodynamic interactions (HI), beads interact with each other and with the walls

94 of the confining geometry. The mobility tensor M, will be discussed in detail in Sec. 6.3.3. For systems with configuration-dependent mobility tensors, the third term in Eq. 6.4 has to be included to obtain the correct dynamics. This is because the configuration dependence of the mobility tensor results in a mean drift with a mean velocity ∂ ∂R

· M.[61] In complex geometries, the expression for M is lacking. Thus to solve Eq.

6.4, we need to turn it into a derivative-free form, and this will also be discussed in Sec. 6.3.3. The tensor B gives the magnitude of the Brownian displacement of the polymer beads, and is coupled to M by the Fluctuation-Dissipation theorem (Eq. 6.5). The vector dW is a random vector composed of independent and identically distributed random variables according to a real-valued Gaussian distribution with mean zero and variance dt. We will use the Chebyshev approximation given in Sec. 6.3.4 to calculate B · dW.

6.3.3

Mobility tensor and Fixman’s midpoint algorithm

To illustrate the components of the mobility tensor, we consider a two-bead chain (a dumbbell) in an unbounded domain, neglecting Brownian effects for the moment. Here the force balance shows that the fluid velocities u(ri ) experienced by bead i are given by 





u∞ 1





1 I ζ

∞

 



G (r1 − r2 )   f1     u1   +  =  , · 1 ∞ G (r − r ) u∞ u2 f I 2 1 2 2 ζ

(6.6)

95 where G∞ (r) =

1 (I 8πηr

+ rr/r 2) is the Oseen-Burgers tensor, the free-space Green’s

function for Stokes’ equations, r = |r|, and I is the identity matrix. The quantity G∞ (r1 − r2 ) · f2 gives the velocity at position r1 generated by the force f2 exerted by the particle at r2 on the fluid. We can write above expression succinctly as, U = U∞ + M · F.

(6.7)

In a confined geometry, a wall correction GW (r1 , r2 ) has to be added to the mobility tensor M, and it becomes 

 M=

1 I ζ

+ GW (r1 , r2 )

∞



G (r1 − r2 ) + GW (r1 , r2 )  . 1 I + G (r , r ) G∞ (r2 − r1 ) + GW (r1 , r2 ) W 1 2 ζ

(6.8)

In general geometries, an analytical expression for GW (r1 , r2 ) is not available, but has to be calculated numerically.[78] In particular, we will see that it is not M that is needed, but rather, the product M · F, which is the velocity generated by the total non-Brownian, non-hydrodynamic forces acting on the beads. Using an accelerated Immersed Boundary Method, we can calculate M · F for polymer in complex geometries efficiently, as discussed in Sec. 6.3.5. Fixman proposed a method to numerically integrate Eq. 6.4 without calculating the derivative term

∂ ∂R

· M.[52] The idea is to approximate the first-order derivative

96 with a finite difference approximation, and transform Eq. 6.4 into R∗ = R +

√

2B(R) · ∆W,

√

2 T ∆R = [U + M(R) · F]∆t + kB T M(R∗ ) · [(B−1 (R)) · ∆W] 2 √ 2 B(R) · ∆W, + 2 ∞

(6.9)

where ∆W is a vector of independent Gaussian random variables with variance ∆t. A derivation of this algorithm is given in Appendix D.

6.3.4

Chebyshev approximation

In Eq. 6.9, there are three matrix-vector multiplications that need to be evaluated: M · x, (B−1 )T · x, and B · x, where x may be any real vector. The Brownian tensor B is coupled to the mobility tensor M by the Fluctuation-Dissipation theorem kB T M = B · BT . A common choice for B is the square root matrix S, satisfying kB T M = S · ST with S = ST . With that, we only need to evaluate M · x and S−1 · x, because (S−1 )T = (ST )−1 = S−1 and S · x = S · (S · S−1 ) · x = kB T M · (S−1 ) · x. Given an algorithm to evaluate M · x, Fixman noted that one can evaluate S−1 · x using Chebyshev polynomial approximation over the range [λmin , λmax ], where λmin and λmax are the minimum and maximum eigenvalues of M respectively [76, 53]. The

97 approximation yL to S−1 · x is a series of matrix-vector multiplications, yL =

L X

al xl ,

(6.10)

0

x0 = dw,

(6.11)

x1 = [da M + db I] · dw,

(6.12)

xl+1 = 2[da M + db I] · xl − xl−1 .

(6.13)

√ Here al are the Chebyshev coefficients of the scalar function 1/ x over the domain [λmin , λmax ], and da = 2/(λmax − λmin )and db = −(λmax + λmin )/(λmax − λmin ). From the above equations, we see that rather than M, all that we need is the product of M and a vector, which can be calculated efficiently by an Immersed Boundary Method as will be described in the next section. To calculate the eigenvalues λ of M, we use power iteration to calculate λmax and Arnold iteration (ARPACK) to calculate λmin .[147] Both algorithms require only dot product calculations between M and some vector x. In summary, with a fast algorithm to evaluate M · x for arbitrary x, we can simulate the Brownian movement of a bead-spring polymer in complex geometries. This fast algorithm is discussed in next section.

6.3.5

Fast Stokes solver with complex boundary conditions

Note that M · F is the fluid velocity generated by the the vector of point forces, F. Thus, to calculate M · F, we need to solve the Stokes flow equations with distributed

98 point forces with no-slip boundary conditions: − ∇p(x) + η∇2 u(x) = − ∇ · u(x) = 0,

X

fn δ(x − xn ),

u(r) = 0, r ∈ ∂Ωb ,

(6.14) (6.15) (6.16)

where p is the pressure, u is the fluid velocity, η is the fluid viscosity, fn is the force exerted on the fluid at point xn , δ(x) is the three-dimensional delta function, and ∂Ωb is the boundary of the fluid domain. In the actual implementation, regularized point forces are used, as further discussed below. Provided a fast Stokes solver for distributed point forces within a periodic domain (GGEM), we developed an accelerated Immersed Boundary Method (IBM) for the Stokes equations with complex boundary conditions by representing the no-slip boundary as forcing term (momentum source) in the Stokes equation. To employ the fast Fourier transform (FFT) technique to speed up the calculation, we chose a periodic domain. We present first the IBM algorithm, then the GGEM algorithm. Immersed boundary method A generic three-dimensional domain is shown schematically in Fig. 6.2. Polymer beads and solid boundaries are embedded in a rectangular periodic domain Ω = [−Lx /2, Lx /2] × [−Ly /2, Lx /2 × [−Lz /2, Lz /2]. The beads of the polymer molecule are represented by the filled symbols, and their position vector is R. The beads exert forces obeying Stokes law on the surrounding fluid. A no-slip boundary ∂Ωb is represented by regularized point forces located at boundary nodes Rb , which are

99

Figure 6.2: Schematic of the immersed boundary method. The system is a three dimensional periodic domain, in this study, a rectangular parallelepiped of Lx × Ly × Lz . The beads of the polymer molecule is represented by the filled symbols, and their position vector is R. The no-slip boundary ∂Ωb is represented by regularized point forces located at boundary nodes Rb , which are uniformly spaced with a boundary mesh size of h.

100 uniformly spaced with a boundary mesh size of h. These boundary nodes do not move relative to the fluid, but rather exert forces on the fluid such that the velocity of fluid at Rb is zero. Hence, the no-slip boundary condition is satisfied at each nodes. The advantages of IBM include its ease of programming due to use of a Cartesian grid, and its potential for parallel processing as the boundary nodes are treated in a similar fashion as the polymer beads. Denoting Mibm as the mobility tensor for all the regularized point forces in the periodic domain, i.e., polymer beads and boundary nodes, the velocities generated by those point forces are, Mibm · Fibm = Uibm .

(6.17)

Partitioning Mibm into smaller tensors for polymer-polymer, polymer-boundary, and boundary-boundary interactions, the relationship between forces and velocities at the polymer and boundary points is given by 









 



pp uni p Mbp   Fp    M  U   U Uibm =  , · +  = b pb bb uni b F M M U U

(6.18)

where p and b denote polymer and boundary respectively, and Uuni is the uniform background velocity which is generated by a constant pressure drop. Here Mpp and Mbb are symmetric and positive definite (and thus invertible) matrices. The quantities Mpb · Fp and Mbp · Fb are the velocities at boundary points generated by polymer forces Fp and the velocities at polymer points generated by boundary forces Fb , respectively, and we are able to calculate them using periodic GGEM, which will be discussed in the next section. Note that the unperturbed Stokes flow velocity U∞ in

101 a complex geometry is the sum of the uniform background flow Uuni and a correcting flow generated by the boundary points Ubr , U∞ = Uuni + Ubr . From Eq. 6.18 we obtain: Up = Uuni + Mpp · Fp + Mbp · Fb ,

(6.19)

Ub = Uuni + Mpb · Fp + Mbb · Fb = 0.

(6.20)

Noting that Mbb is invertible, we obtain: Fb = (Mbb )−1 · (−Uuni − Mpb · Fp ),

(6.21)

Up = Uuni + Mbp (Mbb )−1 · (−Uuni ) + Mpp · Fp − Mbp (Mbb )−1 Mpb · Fp = Uuni + Mbp (Mbb )−1 · (−Uuni ) + [Mpp − Mbp (Mbb )−1 Mpb ] · Fp . (6.22) Comparing the above equation to U∞ = Uuni + Ubr , we obtain Ubr = Mbp (Mbb )−1 · (−Uuni ).

(6.23)

By construction, Eq. 6.20 and Eq. 6.21 also suggest an algorithm to calculate Up : (1) Given the polymer configuration R(tn ), evaluate polymer forces Fp (tn ) including spring forces, excluded volume interactions, and repulsive wall-bead interactions; (2) Calculate the perturbed velocities Upb = Mpb · Fp generated by polymer forces at the boundary points; (3) Determine the boundary forces Fb (tn ) that satisfy Mbb · Fb (tn ) = −Uuni − Upb (tn ) (Eq. 6.21); (4) Evaluate the polymer bead velocities Up = Uuni + Mpp · Fp + Mbp · Fb .

102 Step (3) in this algorithm requires an efficient iterative scheme to calculate the boundary forces that specify the velocity field in the domain we are interested in. In other words, our goal is to solve Mbb · X = U, where U is given and X is unknown. The number of iterations required to reduce the relative error to meet some specified tolerance is a function of the condition number of Mbb (Figure 6.3(a)), which is defined as the ratio of the largest over the smallest eigenvalues of Mbb . We found that a simple conjugate gradient method [123] is sufficient. First, in the algorithm of conjugate gradient method, only matrix-vector multiplications are calculated. Second, it can reduce the relative residual ||Mbb · Fb − U||/||U|| (Eq. 6.21) to less than 10−4 in tens of iterations as shown in Fig. 6.3(b), even when the condition number of Mbb is very large. It should be noted that a similar approach, in a non-Brownian boundary element context, has been proposed by Zhao et al. to study the flow of blood cells in complex geometries.[159] Periodic GGEM The core of our algorithm is a fast solver for the Stokes equations with distributed point forces. In this section, we present the general geometry Ewald-like method (GGEM), which can serve as a “blackbox” to evaluate M · x for a given point force distribution in periodic domain.[67, 68, 66, 114] For periodic boundary condition, Ewald sum and particle-mesh Ewald (PME) methods are available, which are based on the Hasimoto’s solution for Stokes flow driven by a periodic array of point forces. These methods reduce the operation to O(NlogN), compared with O(N 2 ) operations of direct evaluation of all the pairwise interactions.[65, 132, 124] Fast multipole method is another alternative, and the complexity is O(N).[145] GGEM is also an

103

Figure 6.3: Properties of the Mbb tensor for the case where a spherical boundary (R = 3) is represented by uniformly distributed points on the surface and the parameters for the calculation are Lx = Ly = Lz = 10, dx = dy = dz = 0.25, ξ = 4.0. (a) Condition numbers as a function of number of points on the sphere shell. (b) Relative residual as a function of number of iterations.

104 O(NlogN) method when fast Fourier transforms are implemented. The central idea that makes GGEM a fast algorithm is separating the point force density ρ = Σfn δ(x − xn ) into a local part ρl = Σfn (δ(x − xn ) − g(x − xn )) and a global part ρg = Σfn g(x − xn ) by introducing a screening function g(x − xn ). In this work, the screening function g(r) is a modified Gaussian with α as the screening parameter 2 r2 )

g(r) = (α3 /π 3/2 )e(−α

(5/2 − α2 r 2 ),

(6.24)

where r is the radial position vector and r = |r| (Fig. 6.4). This smooth function Z satisfies g(r)dr = 1 for any value of α, and when α → ∞ it becomes a threedimensional Dirac delta function. Correspondingly, the velocity field generated by ρ is

separated to a local contribution and a global contribution u = ul +ug . We chose g(x− xn ) so that an analytical expression for ul is available and ul decays exponentially fast. Thus, the interactions due to ρl are short ranged, and only interactions within some cut-off distance Rc need to be considered in the calculation of ul . In the present case, the global velocity ug is going to be calculated by solving the Stokes equations on a structured mesh (grid) using three-dimensional fast Fourier transform. Note in contrast to conventional Ewald-based methods, FFT is not necessary to achieve computational efficiency in GGEM, so the GGEM approach is not limited to periodic domains.[67, 69, 148, 114] Local velocity field Consider first the local velocity at x, ul (x), which results from the local force density ρl , ul (x) =

X n

Gl (x − xn ) · fn ,

(6.25)

105

Figure 6.4: Screening function g(αr) (Eq. 6.24) for the GGEM algorithm. where the “local Green’s function” is

Gl (x) =

xx erfc(αr) 1 xx 2α 1 2 2 (I + 2 ) − (I − 2 ) 1/2 e(−α r ) . 8πη r r 8πη r π

(6.26)

Here, r = |x| and I is the identity matrix. In this work, to avoid singularity associated with delta point forces, we use regularized point forces by replacing δ(r) with g(r) (Eq. 6.24) with a regularizing parameter ξ and the corresponding regularized local Green’s function is GR l (x) =

1 xx erf(ξr) erf(αr) 1 xx 2ξ 2α 2 2 2 2 (I + 2 )( − )+ (I − 2 )( 1/2 e(−ξ r ) − 1/2 e(−α r ) ). 8πη r r r 8πη r π π (6.27)

When x → 0, 0 GR l =

4α 1 4ξ ( √ − √ )I, 8πη π π

(6.28)

and the velocity at x = 0 stays finite. From Eq. 6.27, we see that because of the presence of the screening function, velocity decays exponentially over a distance pro-

106

Figure 6.5: (Free space) Comparison of the x component of the velocity field driven by a delta point force (lines) acting along +x-axis and by a regularized delta point force (symbols) with a form of modified Gaussian (Eq. 6.24, α = 2.0). (Top) along x-axis. (Bottom) along y-axis. portional to 1/α, provided that α ≪ ξ. Calculation of the local velocity field begins by identifying point forces located within the cut-off distance Rc = 4/α for a given x. Any interactions beyond Rc is ignored. The value of α is empirically determined, and is a constant for certain simulation. To make sure the local velocity satisfies the periodic boundary conditions, we apply the minimum image convention which states that, in the periodic system, the cutoff distance must be smaller than half the shortest box length. In Fig. 6.5, we compare the velocity generated by a delta point force and a regularized point force. The regularizing parameter ξ is different for the polymer beads and the boundary nodes. Consider first the boundary nodes. As with the conventional IBM [110, 97], the boundary force density is distributed onto the grid, which is uniformly spaced with a mesh size of h, through a regularized delta function. In this work, we choose the modified Gaussian as in Eq. 6.24, and we choose ξ such that ξh = O(1). This

107 ensures that the force density associated with each grid node is spread over the length scale of the associated elements, thereby preventing fluid from penetrating the solid boundaries and also preventing unphysically large fluid. For the polymer beads, the choice for ξ is straightforward and is inversely proportional to the hydrodynamic radius of the polymer bead, ξa ≈ 1. In the present work, we choose ξ so that √ ξa = 3/ π, with which the maximum fluid velocity is equal to that of a particle with radius a and the pair mobility remains positive-definite.[69] For any point x in the fluid that is not on the boundary or a polymer bead, the local contribution to the velocity is given by

ul (x) =

n X

GR l (x

i

− xi ) ·

fib

+

m X j

p GR l (x − xj ) · fj .

(6.29)

Here, n and m are the number of boundary points and polymer beads, respectively, which lie within the cut-off distance centered at x. For a boundary point, the local contribution is

ul (xk ) =

n X

GR l (x

i

− xi ) ·

fib

+

m X j

p GR l (x − xj ) · fj .

(6.30)

Note that in the first term, the i = k term requires evaluation of GR l (0), which is given by Eq. 6.28. For a polymer bead, the local contribution is

ul (xk ) =

n X i

GR l (x

− xi ) ·

fib

+

m X

j,j6=k

p GR l (x − xj ) · fj .

(6.31)

The exclusion of the “self-term”, j = k, in the polymer local velocity calculation is

108 because in the Brownian dynamics simulation polymer beads do not move with the fluid velocity. Thus the velocity “seen” by a polymer bead should not include the velocity generated by the bead itself as it moves through the fluid. Global velocity field The velocity field driven by the global part of force density ρg is obtained by solving the Stokes equations with periodic boundary condition on a grid using fast Fourier transform (FFT). The continuous force density g(r) has to be replaced by a grid based force density. In this work, we use the modified Gaussian function to distribute the force density onto the grid points which are within the cut-off distance Rc . This is simply a Fourier collocation approach to the global problem.[27] Again, we use periodic boundary condition with minimum image convention to make sure the forcing term is periodic in all three spatial directions. The linearity of Stokes equation allows us to exploit the efficiency of FFT. We assume the velocity, forcing term, and pressure gradient are periodic in all three spatial directions. The velocity vector U(x) and pressure gradient ∇p(x) at the grid point x = (x1 , x2 , x3 ) are

u(x) =

P

∇p(x) =

ˆ k e−ik·x , u

k

P

k

p ˆk e−k·x .

(6.32)

Here, k = (k1 , k2 , k3 ) is the wavenumber vector and k 2 = k12 + k22 + k32 . Now we apply Fourier transform to the Stokes equations Eq. 6.15 and obtain ˆ k − ηk 2 u ˆ k = −ˆ p f,

ˆ k = 0. ik · u

(6.33)

109 When k 6= 0, taking the scalar product of the Fourier transformed momentum equation with ik yields (k · ˆ f) k, 2 k ˆf (k · ˆf ) − k. = ηk 2 ηk 4

ˆk = − p u ˆk

(6.34)

ˆ 0 = −ˆ When k = 0, p f , which means that the force exerted on the fluid is balanced by the mean pressure gradient of the field. Then the global velocity field is obtained by inverse FFT. To get the velocity at the polymer points from the velocity on the mesh, we perform back-interpolation using quadratic Lagrange polynomials. It is this interpolation step that ultimately controls the order of accuracy of the solution. Finally, note that in the global velocity calculation for polymer beads, we also have to exclude the self-interaction term. This velocity field is determined by the free-space regularized Stokeslet (Eq. 6.26). Therefore, the total velocity seen by a polymer bead is u(xpj ) = ul (xpj ) + ug (xpj ) − lim Gl (x) · fjp . x→0

(6.35)

Validation: GGEM We investigated the accuracy of GGEM as function of screening parameter α and mesh size ∆x, against Hasimoto’s solutions of spatially periodic Stokes equation[65, 113]. The simulation box is a cube of side length L = 10. A delta point force acting along x direction with unit strength is placed at the center of the box. Fig. 6.6 shows Hasimoto’s solution and the numerical solution using GGEM for the x component of velocity at line x = 0, y = 0 and line y = 0, z = 0. Fig. 6.7(a) shows the error as a function of screening parameter α for three different

110

Figure 6.6: Comparison of Hasimoto’s result[65] (symbols) and numerical solution (lines) for the x component of velocity due to point force acting along +x direction in a periodic domain. mesh sizes ∆x. Notice that, empirically, when 20/L < α < 0.8/∆x the error does not depend on α. The lower bound is related to the restriction that cut-off radius Rc = 4/α should be smaller than half of the box size, and the upper bound is about the number of grid points which need to be included to resolve the length scale α−1 . Fig. 6.7(b) shows the error as a function of mesh size ∆x together with least squares fitting to estimate the order of accuracy. The error is defined as the L2 norm of the difference between GGEM results and Hasimoto’s ||E||2 = ||uGGEM − uHasimoto ||2 at 200 sampling points which are randomly distributed within the simulation box. With Lagrange back-interpolation, GGEM has order 3 of accuracy. This result does not change with the positions of the sampling points. Unless otherwise noted, throughout the rest of this study, we choose the screening parameter which satisfies α∆x = 0.8, and ∆x is determined by the desired error tolerance.

111

Figure 6.7: Error ||E||2 as a function of (a) screening parameter α and (b) mesh size ∆x for a point force at the center of a cubic periodic domain.

112 Validation: IBM We now show numerical results of three test cases using the IBM algorithm proposed above: pressure driven flow through a slit geometry (Fig. 6.8), uniform flow around a sphere (Fig. 6.9(a)), and flow generated by a Stokeslet near a solid plane (Fig. 6.9(b)). We use the slit case to investigate the numerical error of IBM. For laminar flow in a slit geometry, the velocity only depends on the perpendicular distance away from the wall, and it is parabolic Ux (z) = 3/2Uavg (1 − (z/H)2 ). In the numerical calculations, the periodic domain is a cube with side of L = 40, and the wall is represented by a square lattice with grid size of h. In all the calculations, we used a mesh size of ∆x = 0.05L for the global velocity calculations. Fig. 6.8(a) shows excellent agreement between the analytic solution and the numerical solution obtained by IBM. The absolute error ||Unumerical − Uexact ||2 is shown as a function of h in Fig. 6.8(b); decays roughly as h1.2 , consistent with the general result that IBM is first-order accurate.[86] The local error is the largest near the boundary and decays very rapidly into the flow (< 1 × 10−5 on the center line). The uniform flow around a sphere case demonstrates that our algorithm works well for curved boundary as well. The exact solution to the flow field in an unbounded domain is: 1 a 3 a Ur = U ∞ [1 − ( ) + ( )3 ]cos(θ), 2 r 2 r 3 a 1 a ∞ Uθ = −U [1 − ( ) − ( )3 ]sin(θ). 4 r 4 r Here, θ is the angle away background flow direction, a = 1 is the radius of the sphere, and r is the radial distance from the center of the sphere. We use a large simulation box (L = 500) to minimize the effect of periodic boundary. When the background flow

113

Figure 6.8: (a) Velocity profile of laminar flow through a slit. The line represents the analytic profile and the symbols represent the numerical result by IBM. (b) Error of IBM as a function of boundary grid size h and the slope is 1.22.

114 is U∞ =(1, 0, 0), we use N = 200 regularized points, which are uniformly distributed on a spherical shell of radius 1, to represent an unit sphere placed at the center of the simulation box. The mesh size is ∆x = 1 and α∆x = 0.1 in this case. From Fig. 6.9(a) we see that IBM gives a satisfactory result. In the case of a Stokeslet above a solid wall, numerical results are tested against the exact results from the analytical expression of the mobility tensor for a point force near a single solid wall, obtained by Blake[22]. The calculation is similar to that in the slit case except that a point force is placed near a wall, and the separation between the walls is large to minimize the effect from the other wall. The parameters for this calculation are L = 40, ∆x = 0.5, and h = 0.5. The wall is the x − y plane, and the Stokeslet is situated at (0, 0, 2), pointing in the x direction. From Fig. 6.9(b) we see that IBM gives excellent results for this case, which involves both force points and boundary points.

6.4

DNA flowing across an array of nanopits

In this section, we apply the immersed boundary method described in last section to study the dynamics of a DNA molecule flowing through a nanoslit with embedded nanopit arrays (Fig. 6.10). The schematic of the problem is given in Fig. 6.10. The device is periodic in the streamwise (x) and spanwise (y) directions. One periodic or one unit cell is a rectangular parallelepiped of Lx × Ly × 2H = 40a × 40a × 8a. Here, a = 77 nm is the bead radius in the bead-spring model representing the DNA. One pit is at the center of the unit cell and has a size of Ld × Ld × H = 16a × 16a × 4a. The pit-to-pit distance is 2Lc = 24a.

115

Figure 6.9: Validation of IBM algorithm. Open symbols are numerical results and lines are analytical calculations. (a) Velocity profile Ux around a unit sphere with no-slip boundary condition. Open circles represent points along the x-axis and open squares are for points along y-axis. (b) Stokeslet near a single wall. Velocity profile of Ux along +z-direction for various (x, y) pair. Open squares are for (x, y) = (3, 3) and open circles are for (x, y) = (1, 1).

116

Figure 6.10: Schematic and the immersed boundary representation of the nanopit problem. Gray spheres indicate points at which no-slip boundary conditions are satisfied. In our simulation, there are two unit cells in a three dimensional periodic simulation box of 2Lx × Ly × 2H. The boundary is represented by boundary points which are regularly spaced with a mesh size of 1a, as shown in Fig. 6.10. Correspondingly, we chose α = 0.8a and ξ = 2.5a as the input parameters for the immersed boundary method. A Peclet number, the ratio between the diffusive and convective time scale, is defined as P e =

Up H , D

and the Weissenberg number is W i = τ0 γ˙ = τ0 UHp , where Up

is the maximum x component of the unperturbed velocity in the y − z plane, D is the chain diffusivity, τ0 is the longest relaxation time of the chain in bulk solution, and γ˙ is shear rate in the pit. Note that P e and W i are linked by the confinement ratio Rg /H, as W i ∼ P e( RHg )2 . For λ−DNA, which has a radius of gyration about 10a, the effective confinement ratio is about 5, considering that the cut-off distance for the excluded bead-wall interaction is 1a. For λ−DNA we use the τ0 value obtained by Jendrejack et al.[74]. Figure 6.11 shows several streamlines of the Stokes flow and the contour plot of the streamwise velocity in the absence of the DNA molecule.

117

Figure 6.11: (Top) Streamlines in the nanopit and (Bottom) contour plot of the streamwise velocity on the x − z plane.

6.4.1

Dynamics at low Peclet number

In the introduction we described the experimental observation of the chain hopping between pits, and the mechanism of a free end diffusing and being dragged downstream from one pit to the next. This behavior is reproduced in our simulations, as shown by the sequence of snapshots in Fig. 6.12. To study the dynamics quantitatively, we tracked the center-of-mass of the DNA molecule and extracted dynamic properties of DNA from its trajectories. Fig. 6.13(a) is a typical time series plot of the x-component of the center-of-mass xc of a λ-DNA molecule. This figure clearly demonstrates that DNA hops from pit to pit at low Peclet number: the molecule stays in one pit for a time interval td , which we will call the residence time, and then jumps to the next pit. As seen in Fig. 6.13(a) at low Peclet number, the time for the

118

Figure 6.12: Snapshots of a hopping event (from (a) to (f)) (top-down view, P e = 3.5).

119 molecule to get to the next pit is very small compared to the residence time. The probability distribution of td is well-fit to an exponential ρ(td ) = 1/τ exp(−td /τ ), as shown in Fig. 6.13(b). This result indicates that, to a good approximation, the hopping events are independent of one another—the hopping dynamics is a Poisson process at low P e, and the rate parameter which characterizes the Poisson process is 1/τ , where τ is the mean residence time. Intuitively, one expects that the larger the applied pressure gradient, the faster the chain hops, or the smaller τ becomes. Fig. 6.14 shows the Peclet-number dependence of the mean residence time τ , nondimensionalized by the longest relaxation time τ0 of λ−DNA. Over the range of Peclet number studied, the results appear to follow a exponential distribution τ ∼ τh exp(−αP e), where 1/τh is basically the hopping frequency as P e → 0. The physical mechanism of hopping becomes clear when we write the Peclet number as P e =

(Up ζc )H , kB T

where ζc is the friction coefficient of the

chain and we have written the chain diffusivity D as

kB T ζc

. The term Up ζc is a measure

of the hydrodynamic drag force on a chain in the pit. Therefore, the product in the numerator is the work done by the fluid to drag the molecule out of the pit, and P e is the ratio of that work to thermal energy. Hence, our simulation demonstrates that at low P e, the dynamics of the hopping is indeed thermally activated. To examine the effect of HI on the dynamics, we also performed simulation for the “free draining” (FD) case, where HI were ignored ( in Fig. 6.14). In FD simulation, the beads are not coupled through hydrodynamical interactions. The mobility tensor is reduced to the product of ζ1 and the identity matrix. Comparing the results obtained for HI and FD, we can see that the mean residence time is indeed affected by the hydrodynamic interactions: it takes a longer time for the molecule to hop to the

120

Figure 6.13: (a) Time series plot of the x-component of the center of mass of a λ− DNA driven through a nanopit array for two low values of Pe. (b) Residence time distribution and best fit to an exponential distribution for conditions shown in (a).

121

Figure 6.14: Mean resident time τ v.s. Peclet number P e. next pit in the FD case. This indicates that although HI is screened on a length scale larger than the slit height, it plays an important role in dynamics happening on smaller length scales. In this case, the “cooperative” motion of DNA segments due to HI in the complex geometry gives rise to an increase of the hopping frequency. We also noticed that the difference between HI and FD gets smaller as the Peclet number increases. As Peclet number increases, the transport process is dominated by DNA-wall steric interaction while hydrodynamic interactions play a relatively minor role. This observation is in line with the observations of Hernandez-Ortiz et al. in the problem of hydrodynamic effects on polymer translocation through a pore.[66] We also studied the relationship between mean residence time and chain length. In general, separation of a mixture of chains of different sizes can only occur if the mean polymer velocity through the channel strongly depends on the chain length. Fig. 6.15 shows the dependence of the mean residence time on chain length for two different

122

Figure 6.15: Mean residence time τ v.s. chain length N and Peclet number P e Peclet numbers. First, as we increase the Peclet number from 3.5 to 5.4, τ becomes smaller. Second, at a given P e, we observe a strong molecular weight dependence of τ at low P e, while τ saturates in the long chain limit. The first observation is consistent with the results shown in Fig. 6.14 that increasing P e lowers the energy barrier height. The mechanism responsible for the second observation is not clear. Based on the barrier crossing theory, we have the following hypothesis. For a given pit size, there is a corresponding DNA size Lc , which fills the pit. For a DNA molecule longer than Lc , the energy barrier is determined only by the chain segments with the same length as Lc , i.e. it is determined by the part of the chain residing in the pit. Since the transport velocity of DNA in the slit is independent of chain length [141], we can assume the reaction rate prefactor in the barrier crossing theory is also independent of chain length, thus, in the long chain limit, the mean residence time saturates. The effect of HI gets more pronounced as the chain length gets longer, as shown by the

123 empty square symbols in Fig. 6.15. This is consistent with the observation of Izmitli et al.[73] and Hernandez-Ortiz et al.[66] in the problem of hydrodynamic effects on polymer translocation through a pore.

6.4.2

Dynamics at high Peclet number

When we increase the pressure drop, the transport characteristics of DNA in the device change, as indicated by the change of the shape of the probability density function of residence time (Fig. 6.16(a)). In contrast with ρ(td ) at low Peclet number, the distribution at high Peclet number can be fitted to a Gaussian ρ(td ) = 1 √ w 2π

2

−Tc ) ), with mean Tc and variance w 2 . This shift from exponential to exp(− (td2w 2

Gaussian reflects the change of the dominant physics as P e increases, from a primarily stochastic process at low Peclet number to a primarily deterministic one at high Peclet number. Qualitatively, our simulation results are consistent with experimental observations [36]. We do not seek a quantitative comparison, as data about the mean residence time were not reported in the experimental study. Fig. 6.16(b) shows the mean residence time as function of P e, which is also quite different from that at low P e. It appears to follow a power law at high Pe, Tc ∼ 1/P e. This is not surprising as at high Pe, as the deterministic drifting process dominates. Hence the residence time is well-approximated by the time it takes to travel across one pit, which is approximately Tc ≈ Lx /Up ∼ 1/P e. The variance of the residence time, w 2 , approaches a scaling of P e−2 in the high Peclet number limit. This is consistent with the macrotransport theory prediction of force-driven transport of a point-size Brownian particle in a slowly varying periodic channel.[25, 84]

124

Figure 6.16: (a) Residence time distribution at high Peclet number and best fit to a Gaussian distribution. (b) Fitting parameters for Gaussian distribution at different Peclet numbers.

125

6.5

Conclusions

In this work, we presented an immersed boundary method (IBM) based on the general geometry Ewald-like method (GGEM), aimed at simulating Brownian motion of polymer in complex geometries with full description of hydrodynamic interactions. The core idea is to replace the complex boundary with regularized point forces, chosen to satisfy the appropriate boundary conditions. The IBM-GGEM methodology correctly reproduces the velocity field in several test cases including cases involving curved boundary and point momentum sources. As an application of this methodology we have considered the simulation of a single DNA molecule driven through a nanoslit with embedded nanopits array, by a pressure gradient. We studied the dynamics of the DNA as a function of the Peclet number and chain length, as well as the influence of hydrodynamic interactions by comparing with free draining simulation results. We found that the transport characteristics of the hopping dynamics in this device differ at low and high Peclet number, and for long DNA, relative to the pit size, the dynamics is governed by the segments residing in the pit. We also found that HI plays an important quantitative role in the hopping dynamics even in such a highly confined system. We expect that our algorithms will find many applications in micro- and nanofluidics.

126

Chapter 7 Conclusion and future work In conclusion, we have developed general numerical frameworks to study electrophoreticdriven and flow-driven DNA dynamics in complex geometries, and we presented a systematic investigation of three different problems using Brownian dynamics simulations. In both the problem of soft nanomechanical elements and the problem of flowing DNA through nanopits, the concepts of entropic trapping and barrier crossing give satisfying explanations to the observed dynamics. It seems that we have a good understanding of the basic physics in those microfluidic devices, where DNA dynamics is determined by the combined effects of major driving forces and chain configuration entropy. Now the goal is shifted to apply the known physics to the design of devices. This echoes the motivation of this work, which is to build validated numerical tools to help the design and optimization of novel processes and devices in a very broad area. Considering the future of both the simulation methods and their applications

127 to the problem of DNA dynamics in complex geometries, harnessing the increasing computational power to conduct more detailed simulations of larger system for longer time is an inevitable trend. Future work may include: • New mechanical model for DNA under nanometer scale confinement. When a DNA chain is confined to a space smaller than its persistence length, which is already achieved in experiments, the bead-spring model breaks down. We thus need to build a finer-scale model of the DNA to study these systems, for example, a bead-rod model with a rod length of 10 nm. For a λ-DNA molecule, this means thousands of beads, which is out of the reach of traditional Brownian dynamics simulation technique when long-ranged hydrodynamic interactions are included. But it might be feasible with our new accelerated immersed boundary method using graphic processing unit. • Graphic processing unit (GPU) acceleration. The parallel architecture of GPU makes it very effective for algorithms where processing of large blocks of data is done in parallel. It has been applied to accelerate molecular dynamics simulation and computational fluid dynamics calculation. As recently reported, GPU calculations achieved tens of fold, even hundreds of fold in some case, speedup over an optimized CPU implementation. The advantages of our accelerated Brownian dynamics algorithm is that it divides the long-ranged interactions into short-ranged pair-wise interactions, which can be calculated using analytical expressions, and long-ranged interactions, which change smoothly with small gradient and can be calculated on a regular structured Cartesian grid. Both calculations are straightforward to parallelize efficiently. Also, in our algorithm,

128 the boundary nodes are treated in a similar fashion as the polymer beads, which gives its potential for parallel processing.

129

Appendix A Lattice random walk model of a tethered polymer In this appendix, we use a random walk lattice model to calculate the thermodynamic properties of an ideal polymer molecule tethered to a hard nonadsorbing wall. For an ideal chain, the entropy S(r) associated with all chain conformations starting from an origin and ending at r is simply related to the number of distinct walks Z(r) connecting the origin and r in n steps.[35]

S(r) = kB ln[Zn (r)]

(A.1)

Therefore, the free energy of an ideal chain is F (r)/kB T = − ln[Zn (r)]. It is useful to know how does the free energy change as a function of a single reaction coordinate, such as the perpendicular distance of the free end to the wall ze for a tethered molecule. It could reduce the problem to one dimensional and allow us to apply Kramer’s theory of barrier crossing to study the dynamic properties of the walk. Here we introduce

130 the potential of mean force of a polymer chain as a function of ze F (ze ) = F ∗ (ze ) − kB T ln[

hρ(ze )i ] hρ∗ (ze )i

(A.2)

where F ∗ (ze ) and hρ∗ (ze )i are arbitrary reference values.[121] The average distribution function hρ(ze )i along the reaction coordinate is obtained from a Boltzmann weighted average hρ(ze )i =

R

′

drδ(ze (r) − ze )e−U (r)/kB T R , dre−U (r)/kB T

(A.3)

where U(R) represents the total energy of the chain as a function of the beads’ ′

positions R and δ(ze (R) − ze ) picks out the configurations with the same ze . For an ideal chain, every path has the same weight, and above expression reduces to

hρ(ze )i =

Zn (ze ) , Zntot

(A.4)

and therefor, the potential of mean force is F (ze ) = F ∗ (ze ) − kB T ln[

Zn (ze ) ]. Zn∗ (ze )

(A.5)

Above expressions indicate one can calculate the free energy profile along the chosen reaction coordinate (ze ), by calculate the ratio of the number of configurations of a specific state to that of the reference state. The property of an ideal chain is often studied using a lattice random walk model [35]. For a n-step random walk on an integer lattice, we can directly count the number of walks starting origin and ending at (x, y, z). When the walk is near a wall, absorbing boundary condition should be used, i.e., we should discard those walks

131 touching the wall at least once.[40, 134] Once the number of walks satisfying both absorbing boundary and various constraints are obtained, the free energy is readily obtainable using thermodynamic relationships presented above. Next, we present the results of combinatorics calculations for counting the chain configurations. According to Bertrand’s “ballot theorem”, the number of 1-D random walks on the integers of n steps, from the origin to the point x > 0 and never return to the origin, is Zt1d (x, N)

x = N



 N +x , N

(A.6)

2

assuming N and x have the same parity mod(N − x, 2) = 0. Thus, the total number of walks is Zt1d

=

X

Wt1d (x, N)

x

X x N  = N +x . N 2 x

(A.7)

Assuming N = 2k + 1(k = 0, 1, ...), Zt1d

  X x N  1 X N = = x N +x N +x N N 2 2 x=1       x=1  1 N N N N = [ ]. + (2k + 1) + ... + (2k − 1) +3 k+k+1 k+k k+2 N k+1 (A.8)

Note that (2i − 1) = (k + i) − (k − i − 1),

N k+i




=

N k+1−i




and (k + i)

N k+i




=

132 (N − k + 1) Zt1d

N N −k+1




. Last expression is

        1 N N N N = [ ]. + (2k + 1) + ... + (2k − 1) +3 k+k+1 k+k k+2 N k+1         N N N N 1 + N] − + ... + (2k) −k = [(k + 1) 1 2k k N k+1     1 2k N = (k + 1) . (A.9) = k k+1 N

Similarly, we can prove that when N = 2k(k = 1, 2...), Zt1d =
N −1 Zt1d = [N/2] , where [x] is the largest integer smaller than x.

2k−1 k


. Therefore,

In 2-D case, assuming there are i steps in the wall normal direction and N − i

steps in the wall parallel direction, the total number of paths are

Zt2d

    N  X 2N − 1 N − 1 i − 1 N −i . 2 = = i N − 1 ] [ i − 1 2 i=1

(A.10)

Similarly, in 3-D, the total number of paths is

Zt3d

  N  X N − 1 i − 1 N −i = 4 . i [ i − 1 ] 2 i=1

(A.11)

In summary, we obtain the analytical expressions of the number of configurations of

133

x

z

y N steps 2 2

1

0

Figure A.1: A tethered polymer as lattice random walk. a tethered polymer molecule, and the total numbers of paths are     N X N −1 2N x N √ , = ≈ = N N +x [ ] N 2πN 2 2 x=1     N  X N − 1 m − 1 N −m 4N 2N − 1 2d √ Zt (N) = , ≈ 2 = m [ m − 1 N − 1 ] 4πN 2 m=1   N  X 6N N − 1 m − 1 N −m 3d √ , 4 ≈ Zt (N) = m ] [ m − 1 6πN 2 m=1 Zt1d (N)

(A.12) (A.13) (A.14)

where m is the number of steps in the confined(wall normal) direction. There is √ no closed expression for the sum in 3-D. However, the approximation 6N / 6πN converges to the exact result as N increases (data not shown), following the same scaling as those for 1-D and 2-D cases. In this work, the total number of N-step tail is called the fundamental solution of a tethered chain, and it is denote Z0N . When a polymer is grafted to a bilayer, it can induce gelation or other phase

134 changes in lamellar phases, which might be related to the interactions of locally deformed membrane patches which are created by the entropic forces exerted by the grafted polymer molecules[20]. When an ideal polymer chain is grafted to a hard wall, the number of accessible configuration is greatly reduced which leads to an entropy penalty. Phenomenologically, we observe a repulsive force which drives the polymer away from the wall, i.e., we need to apply a tethering force to fix the end of polymer. We are interested in the dependence of this tethering force on polymer chain length. For a lattice polymer, the required force f0 is related to the number of configurations Z as f0 = −

F0 − F1 kB T ∂F (l) |l=0 ≈ = (ln Z0N +1 − ln Z0N ), ∂l 1a a

(A.15)

where l is the perpendicular distance of the tethered end from the wall, the lattice grid size is a = 1 which has a physical meaning of persistence length of the polymer molecule. Using the expression of the number of tail configuration (Eq. A.14), we know that Z0N ≈

N √6 6πN

N+1

and Z0N +1 ≈ √ 6

6π(N +1)

. Plug these into last equation and we

obtain kB T Z0N +1 ln N a Z r0 N kB T 1 N +1 kB T ln 6 = (ln 6 − ln ) ≈ a N +1 a 2 N kB T 1 ≈ (ln 6 − ). a 2N

f0 =

(A.16)

Clearly, the force is repulsive, and the longer the chain, the larger the force. Interestingly, the force becomes independent of the chain length in the long chain limit. The expression we obtained is similar to that of Bickel et al. [20] using the

135 1

10

0

∆F(l/Rg)/kBT

10

−1

10

−2

10

N=6 N=8

−3

10

−4

10

N=10 N=12

−5

10

0

1

2

l/Rg

3

Figure A.2: The wall-polymer potential calculated using the fundamental solution of a tethered chain for various chain length (symbols). Dashed line is a fit to Eq. A.23 diffusion equation[20], f0 =

kB T a

2

a exp(− 4R 2) ≈ g

kB T a

(1 −

a2 ), 4R2g

and it has the same

dependence on persistence length a and on chain length n. The difference of the constant ln 6 is due to the nature of the lattice. To our best knowledge, there is no experimental results for us to compare with. For a semi-flexible fluorescently stained λ-DNA with a = 50nm and Rg = 770 nm, the tethering force is about 0.82 pN. For a flexible molecule of poly(methacrylic acid) (PMAA) with a=0.3nm and Rg = 700nm, however, the tethering force is 11.2 pN. Above analysis can be extended to calculate the effective repulsive potential between a free polymer molecule and a hard wall using the fundamental solution of a tethered molecule. The total number of N-step paths in the half space starting l away

136 from the wall is denoted WlN , and it can be obtained recursively as, W0N = Z0N ,

(A.17)

W1N = Z0N +1 = 4W1N −1 + W2N −1 ,

(A.18)

W2N = Z0N +12 − 4Z0N +1 = 4W2N −1 + W1N −1 + W3N −1

(A.19)

W3N = Z0N +2 − 8Z0N +1 + 15Z0N +1 = 4W3N −1 + W2N −1 + W4N −1

(A.20)

N N Wi+1 = WiN +1 − 4WiN + Wi−1 .

(A.21)

We can see that WlN is a linear combination of the fundamental solutions Z0N +i , i = 1, 2, ..., l − 1, and the coefficients can be constructed recursively. Once all WlN are calculated, taking the free chain as the reference state, the wall-polymer potential can be calculated according to Eq. A.5 as WlN ∆F (l, N) = − ln N . kB T 6

(A.22)

Results for various chain lengths are plotted in Fig. A.2. First, when l is normal√ ized by the size or the molecule Rg = Na, where a is the lattice size, all results collapse onto a master curve. Second, the interaction potential decreases exponentially fast when l < Rg and it can be fitted to the following equation ∆F id (l, N) l = − ln(erf( )), kB T Rg

(A.23)

which is obtained from a continuous description of the random walk using the diffusion equation[24]. This repulsive potential can be used in mesoscale simulation of a polymer molecule near a single wall.

137

Appendix B BD/FEM algorithm for simulating DNA electrophoresis The objective of this appendix is to present a numerical framework, Brownian dynamics (BD)/Finite element method (FEM), to simulate DNA electrophoresis in complex geometries. The proposed numerical scheme is composed of three parts: • a bead-spring DNA model, • a finite element scheme for the non-homogeneous electric field, • and a coupling algorithm to find the values of the electric field at the DNA beads location. The bead-spring model has been described in Chapter 3. In the reminder of this appendix, we discuss the finite element scheme and the coupling algorithm. For a linear DNA molecule, the Debye length (κ−1 ), which is the length scale over the charges on the DNA backbone are screened out by the counter-ions, is much

138 smaller than the persistence length of DNA in sufficiently concentrated salt solution (typically κ−1 is approximately 2 nm in thickness under physiological conditions). Thus, DNA globally behaves like a neutral polymer and the electric field is governed by Laplaces equation, since local electric field disturbance due to a DNA molecule is screened over the Debye length. Also, we consider the system where electroosmotic flow is eliminated using a polymer layer grafted on the walls. In terms of geometry, we consider a microfluidic device with a constant channel height, and the channel height is usually much smaller than the other length scales of the device. Therefore, the problem is reduced to two dimensional. Still, with increasing complexity of geometry, a numerical method is required to calculate the electric field in a microfluidics with complex geometries. The electric field is computed from the electric potential with the relationship E = −∇φ in an arbitrary domain, where φ denotes electric potential. The governing equations for the electrostatic potential are,

φ = φ0

or

∇φ2 = 0 in Ω,

(B.1)

∂φ = 0 on ∂Ω, ∂n

(B.2)

where n is the normal vector at the boundary. We consider two types of boundary conditions. One is Dirichlet boundary condition which specifies the values of potential on the boundary of the domain, such as the electric potentials at the inlet and outlet. The other is Neumann boundary condition which specifies the values that the derivative of potential is to take at the insulating channel walls. In experiments, the device is usually made of poly(dimethylsiloxane) (PDMS) which is an insulator.

139 Hence, a zero gradient boundary condition is assigned along the channel walls. In this work, we use a commercial partial differential equation (PDE) solver, COMSOL, which is based on the finite element method (FEM), to solve the governing Laplace’s equation for the electrostatic potential (Eq. B.2). For detailed information on FEM, the readers are referred to [17] and COMSOL manual. In a microfluidic device, steep field gradient exists around obstacles or sharp corners, which means that the space around an immersed structure like a post, and around a sharp corner like near a contraction/expansion region, should be more finely discretized. Once electric potential is obtained, E is simply computed from E = ∇φ and stored for lookup. To calculate the field strength at the bead location, we need a fast algorithm to find the triangular element surrounding the bead. For an irregular triangular mesh, finding the electric field at a specified point in the domain is not a trivial problem. As we just discussed above, to resolve the steep field gradient around the obstacles and corners, very fine mesh needs to be used. This means that we have a large number of elements to look through. We address this task as follows. We first assume that the target point is closer to the centroid of the element enclosing it, than to the centroids of all the other elements. Of course this is not always true as shown in Fig. B.1. So we need to find the four nearest centroids to the target, and then find the one enclosing it using a confirming algorithm. Now the problem is reduced to find the 4-nearest centroids to the target. We can solve it using the nearest neighbor search algorithm, which is a well-known algorithm in computational geometry [133]. In this work, we use a variant of this algorithm called Approximate Nearest Neighbor Searching (ANN) developed by David M. Mount and Sunil Arya [11]. Finding the nearest point needs O(log N) operation in the case of

140

Figure B.1: Schematic of the nearest neighbor search algorithm (NNS). To find the element enclosing the target point, we use NNS to find the k-nearest centroids of the elements. randomly distributed points. The worst case search time, when a 2-dimensional KDtree containing N nodes is used, is O(N 1/2 ). This allows us to use very fine mesh in the calculation.

141

Appendix C Dimer in shear flow In this appendix we analytically and numerically consider the case of a Brownian particle attached to the origin with a harmonic spring with stiffness k = κζ, where ζ is the friction coefficient, subject to shear flow, vx = γy. ˙ If additionally we include a hard wall at y = 0 such that the random walk is restricted to the upper half-space, y > 0, we are essentially considering a tethered polymer chain composed of two beads. Note that this problem is essentially two-dimensional. The overdamped Langevin equations for the particle coordinates are

x˙ = γy ˙ − κx + Fx y˙ = −κy + Fy , where F denotes the random forcing. After performing a Fourier transform in time,

142 we get the solution in Fourier space (iν + κ)Fˆx + γ˙ Fˆy (iν + κ)2 Fˆy , yˆ = (iν + κ)

xˆ =

D ⋆ E from which we can obtain all cross-correlation functions using the identities Fˆx Fˆx = D ⋆ E D ⋆ E Fˆy Fˆy = α and Fˆx Fˆy = 0. In particular, we obtain the monotonicallydecreasing non-normalized PSD



ˆ ˜ Sxy (ν) = Cxy = khˆ x⋆ yˆik =

(ν 2

αγ˙ , + κ2 )3/2

and, after an inverse Fourier transform, the non-normalized CCF    αγe ˙ −κt (2κt + 1)/(4κ2) for t ≥ 0 ˜ Cxy (t) = .   αγe ˙ κt /(4κ2 ) for t < 0

In the case of no shear flow, γ˙ = 0, we obtain that C˜xx (t) = αe−κ|t| /(2κ), showing that the relaxation time is τ = κ−1 and thus Wi = γ/κ. ˙ For the harmonic spring dimer the relaxation time does not depend on Wi. The cross-correlation function shows a single peak at tmax = (2κ)−1 = τ /2, and after proper normalization, Cxy (t) = q ˜ Cxy (t)/ C˜xx (t = 0)C˜yy (t = 0), the height of the peak in the CCF is found to be max Cxy

√ −1/2 2e Wi . = Cxy (τ /2) = p 2 + Wi2

(C.1)

This analytically-solvable dimer model, even without a hard wall, reproduces the

143 characteristics of the CCF that we observe for tethered polymer chains in shear flow. Specifically, Cxy (t) has an asymmetric peak of width ∼ τ centered at t = τ /2 and height ∼ Wi. There is no periodicity in the motion of the dimer and no “cycling” time-scale other than the intrinsic relaxation time τ . The dimer problem can no longer be solved analytically if a hard wall is present or if the spring is non-linear (e.g., FENE or worm-like). We can, however, study the dimer with a non-linear spring and/or in the presence of a hard wall numerically using Brownian Dynamics (without hydrodynamics). Some results for Wi = 2 are given in Fig. C.1, where we also show the analytical solution for the harmonic dimer and the results for longer tethered chains. When a hard wall is present, the numerical results show that the position of the peak in the CCF shifts to smaller times and reduces in height. For the non-linear springs, the position of the peak moves to smaller times as Wi increases, exactly as we observe for the tethered chains. The height of the peak is several times larger for a dimer than for a chain with N ≫ 1 beads, which is not unexpected. Even after including non-linearity and the hard wall, the dimer model fails to reproduce the smaller but still substantial negative peak at t < 0 that we observe in the CCFs for the longer tethered chains at small Wi. An analytical calculation for a harmonic chain tethered to a point and subjected to shear flow might reproduce that feature as well. We can mimic such a peak by constructing an artificial CCF, Cˆxy (t) = Cxy (t) − αCxy (−t),

(C.2)

where 0 < α < 1 controls the depth of the negative peak, and Cxy is the analytical

144 CCF for the harmonic dimer . As illustrated in Fig. C.1, such an empirical fit matches the numerical results quite well. The Fourier transform of Eq. (C.2) gives an empirical PSD of the form ˆ S(2πν = Ω/τ ) ∼ Wi

p

(1 + α2 )(1 + Ω2 ) − 2α(1 − Ω2 ) , 1 + 2Ω2 + Ω4

which for α > 1/3 exhibits a wide maximum at frequencies Ω = 2πτ /T ∼ 0.5, i.e., at a period T ∼ 10τ . As illustrated in Fig. C.1, the maximum in this PSD is very reminiscent of the “peaks” in the PSD observed in Refs. [127, 37, 38], where they were attributed to the existence of a periodic motion with period of about 10τ . The analytical shape of the PSD only involves τ as a relevant timescale, and the crosscorrelation function has an exponential decay at large times ∼ exp(−t/τ ), just like the autocorrelation function for the end-to-end vector used to define relaxation times. Such an exponential decay is inconsistent with periodic motion, but is consistent with some recent theoretical models that suggest similar correlations for a free chain in shear flow [34, 156]. In summary, as seen from this simple analytical example of a dimer in a flow, a maximum in the PSD does not imply any periodic motion and the claim of an existence of a new physical timescale other than the internal relaxation time of the polymer is not justified.

145

0.40.4 α=1/2 α=2/3 α=3/4

11

0.30.3

PSD PSD

C xy

C xy

0.20.2 0.10.1

0.10.1

0.0 0

Rescaled theory C xy(2τ)/2 Harmonic dumbbell FENE dumbbell Wormlike dumbbell Rescaled chain 3C xz(τ) (BD, N=20)

-0.1 -0.1

Empirical fit ( α=2/3)

-0.2 -0.2 -2 -2

-1

-1

0

0 t/τ

t/t

1

1

22

33

0.01 0.01 0.01 0.01

0.1 Dimensionless frequency f=ντ

0.1

11

Dimensionless frequency f=nt

Figure C.1: The left panel shows the cross-correlation function for a dimer (dumbbell) tethered to a hard wall and subjected to shear flow, for a harmonic, FENE and a worm-like spring. We also show a rescaled form of the analytical solution for a harmonic dimer in shear flow (without a hard wall). The height of the peak diminishes by a factor of about 2 when a hard wall is present, so we have rescaled the analytical solution for the harmonic dumbbell accordingly. The position of the peak shifts to smaller times when a hard wall is present as well, and we have thus rescaled the time for the analytical solution. The CCF for a wormlike chain of N = 20 beads, as obtained from Brownian Dynamics simulations, is also shown for qualitative comparison after scaling by a factor of 3 to bring its height in agreement with the dimer case. We also show an empirical fit to the Brownian Dynamics simulations of the form proposed in Eq. (C.2), for which the PSD can be analytically calculated and shows a maximum at period T ∼ 10τ , depending on the value of the tunable parameter α, as illustrated in the right panel.

146

Appendix D Fixman’s midpoint algorithm Fixman[52] proposed a method that avoids computing the derivative of diffusion tensor, which is generally unknown for complex geometries. As a simple illustration of how this method works, we start with the one dimensional case with only the drift caused by the position-dependence of D, and a forward Euler step is simply

∆x = x(t + ∆t) − x(t) = (

√ d D(x))∆t + 2B(x(t))∆w. dx

(D.1)

We can approximate the derivative using a finite difference, in which case this equation becomes ∆x =

√ D(x(t) + ∆x) − D(x(t)) ∆t + 2B(x(t))∆w. ∆x

(D.2)

This is an implicit method as information at t + ∆t is needed. Fixman proposed an intermediate step defined as ∆x∗ = x∗ (t + ∆t) − x(t) =

√

2B(x(t))∆w.

(D.3)

147 This is a valid approximation, as in Eq. D.2 the first term is proportional to ∆t and the second term is proportional to ∆w ∝ ∆t1/2 . When ∆t → 0, the first term is vanishingly small compared to the second term. Then the finite difference approximation of the derivative term becomes D(x∗ ) − D(x(t)) D(x∗ ) − D(x(t)) D(x(t) + ∆x) − D(x(t)) √ ≈ = . ∆x ∆x∗ 2B(x(t))∆w

(D.4)

Inserting this into Eq. D.2 and notice that ∆t = (∆w)2 when ∆t → 0, we finally have ∆x =

√ D(x∗ ) − D(x(t)) √ ∆w + 2B(x(t))∆w. 2B(x(t))

(D.5)

The vector version of this, with D = B · B, is given by √ 2/2(D(r∗) − D(r(t)))(B−1 (r(t)))T · ∆w + 2B(r(t)) · ∆w √ √ = 2/2D(r∗)(B−1 (r(t)))T · ∆w + 2/2B(r(t)) · ∆w.

∆r =

√

(D.6) (D.7)

We now show that the first term on the RHS of Eq. D.6 equals (∂/∂r · D)dt when √ ∆t → 0. The ith component of the first term is 2/2(Dij (r∗ ) − Dij (r))(B −1 )Tjk ∆wk . Expand Dij around r = r(t) Dij (r∗ ) = Dij (r) +

∂Dij ∗ ∆rl + O(∆rl∗2 )). ∂rl

(D.8)

Plug in Eq. D.3 and ignore the higher order terms: Dij (r∗ ) − Dij (r) =

∂Dij √ 2Blm ∆wm . ∂rl

(D.9)

148 Plug this into √

2/2(Dij (r∗ ) − Dij (r))(B −1 )Tjk ∆wk =

√

2/2

∂Dij √ 2Blm ∆wm (B −1 )Tjk ∆wk ∂rl

∂Dij Blm (B −1 )Tjk ∆tδmk ∂rl ∂Dij −1 = Blk Bkj ∆t ∂rl ∂Dij ∆tδlj = ∂rl ∂Dij ∆t. = ∂rj =

In the second step, we use the property of Wiener process dwi dwj = δij dt.

(D.10)

149

Bibliography [1] M. D. Graham. Fluid Dynamics of Dissolved Polymer Molecules in Confined Geometries. Annu. Rev. Fluid Mech., 43:273, 2011. [2] J. J. de Pablo. Coarse-grained simulations of macromolecules: From DNA to nanocomposites. Annu. Rev. Phys. Chem., 62:555, 2011. [3] R. Adhikari, K. Stratford, M. E. Cates, and A. J. Wagner. Fluctuating Lattice Boltzmann. Europhys. Lett., 71:473, 2005. [4] U. S. Agarwal, A. Dutta, and R. A. Mashelkar. Effect of flexibility on the shear-induced migration of short-chain polymers in parabolic channel flow. J. Fluid Mech., 557:297, 2006. [5] P. Ahlrichs and B. Duenweg. Simulation of a single polymer chain in solution by combining lattice Boltzmann and molecular dynamics. J. Chem. Phys., 111:8225, 1999. [6] B. J. Alder and T. E. Wainwright. Studies in molecular dynamics. I. General method. J. Chem. Phys., 31:459, 1959.

150 [7] F. J. Alexander and A. L. Garcia. The Direct Simulation Monte Carlo Method. Comput. Phys., 11:588, 1997. [8] I. Ali and J. M. Yeomans. Polymer translocation: The effect of backflow. J. Chem. Phys., 123:234903, 2005. [9] A. Alvarez and R. Soto. Dynamics of a suspension confined in a thin cell. Phys. Fluids, 17:093103, 2005. [10] P. E Arratia, C. C Thomas, J Diorio, and J. P Gollub. Elastic instabilities of polymer solutions in cross-channel flow. Phys. Rev. Lett., 96:4, 2006. [11] S. Arya and D. M. Mount. Approximate nearest neighbor queries in fixed dimensions. In Proceedings of the fourth annual ACM-SIAM Symposium on Discrete algorithms, SODA ’93, pages 271, Philadelphia, PA, USA, 1993. Society for Industrial and Applied Mathematics. [12] P. J. Atzberger, P. R. Kramer, and C. S. Peskin. A stochastic immersed boundary method for fluid-structure dynamics at microscopic length scales. J. Comput. Phys, 224:1255, 2007. [13] C. Aust, M. Kroger, and S. Hess. Structure and dynamics of dilute polymer solutions under shear flow via nonequilibrium molecular dynamics. Macromolecules, 32:5660, 1999. [14] A. Balducci, P. Mao, J. Han, and P. S. Doyle. Double-stranded DNA diffusion in slitlike nanochannels. Macromolecules, 39:6273, 2006.

151 [15] J. R. Banavar, M. Cieplak, T. X. Hoang, and A. Maritan. First-principles design of nanomachines. Proc. Natl. Acad. Sci. U. S. A., 106:6900, 2009. [16] A. J. Banchio and J. F. Brady. Accelerated Stokesian dynamics: Brownian motion. J. Chem. Phys, 118:10323, 2003. [17] K. J. Bathe. Finite Element Procedures. Prentice-Hall, Englewood Cliffs, NJ, 1996. [18] V. A. Beck and E. S. G. Shaqfeh. Ergodicity breaking and conformational hysteresis in the dynamics of a polymer tethered at a surface stagnation point. J. Chem. Phys., 124:094902, 2006. [19] J. B. Bell, A. Garcia, and S. A. Williams. Numerical Methods for the Stochastic Landau-Lifshitz Navier-Stokes Equations. Phys. Rev. E, 76:016708, 2007. [20] T. Bickel, C. Jeppesen, and C. M. Marques. Local entropic effects of polymers grafted to soft interfaces. Eur. Phys. J. E, 4:33, 2001. [21] R. B. Bird, C. F. Curtiss, R. C. Amstrong, and O. Hasager. Dynamics of Polymer Liquids Vol. 2. Wiley-InterScience; 2 edition, 2001. [22] J. R. Blake. Image system for a Stokeslet in a no-slip boundary. Proc. Camb. Philos. S-M, 70:303, 1971. [23] B. M. Boghosian, P. J. Love, P. V. Coveney, I. V. Karlin, S. Succi, and J. Yepez. Galilean-invariant lattice-Boltzmann models with H theorem. Phys. Rev. E, 68:025103, 2003.

152 [24] P. G. Bolhuis, A. A. Louis, J. P. Hansen, and E. J. Meijer. Accurate effective pair potentials for polymer solutions. J. Chem. Phys., 114:4296, 2001. [25] H. Brenner and D. A. Edwards.

Macrotransport Processes.

Butterworth-

Heinemann, 1993. [26] C. Bustamante, J. F. Marko, E. D. Siggia, and S. Smith. Entropic elasticity of lambda-phage DNA. Science, 265:1599, 1994. [27] C. G. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang. Spectral Methods: Fundamentals in Single Domains. Springer-Verlag, 2006. [28] C. Chang, C. Liu, and C. Lin. Boundary conditions for lattice Boltzmann simulations with complex geometry flows. Comput. Math. Appl., 58:940, 2009. [29] R. Chelakkot, R. G. Winkler, and G. Gompper. Migration of semiflexible polymers in microcapillary flow. EPL, 91:14001, 2010. [30] Y. L. Chen, M. D. Graham, J. J. de Pablo, K. Jo, and D. C. Schwartz. DNA molecules in microfluidic oscillatory flow. Macromolecules, 38:6680, 2005. [31] Y. L. Chen, M. D. Graham, J. J. de Pablo, G. C. Randall, M. Gupta, and P. S. Doyle. Conformation and dynamics of single DNA molecules in parallel-plate slit microchannels. Phys. Rev. E, 70:060901(R), 2004. [32] Y. L. Chen, H. Ma, M. D. Graham, and J. J. de Pablo. Modeling DNA in confinement: A comparison between the brownian dynamics and lattice Boltzmann method. Macromolecules, 40:5978, 2007.

153 [33] B. Chun and A. J. C. Ladd. Interpolated boundary condition for lattice Boltzmann simulations of flows in narrow gaps. Phys. Rev. E, 75:066705, 2007. [34] D. Das and S. Sabhapandit. Accurate Statistics of a Flexible Polymer Chain in Shear Flow. Phys. Rev. Lett., 101:188301, 2008. [35] P. G. de Gennes. Scaling Concepts in Polymer Physics. Cornell University Press, 1979. [36] J. T. Del Bonis-O’Donnell, W. Reisner, and D. Stein. Pressure-driven DNA transport across an artificial nanotopography. New J. Phys., 11:075032, 2009. [37] R. Delgado-Buscalioni. Cyclic motion of a grafted polymer under shear flow. Phys. Rev. Lett, 96:088303, 2006. [38] R. Delgado-Buscalioni. Dynamics of a Single Tethered Polymer under Shear Flow. AIP Conference Proceedings, 913:114, 2007. [39] R. Delgado-Buscalioni and G. De Fabritiis. Embedding molecular dynamics within fluctuating hydrodynamics in multiscale simulations of liquids. Phys. Rev. E, 76:036709, 2007. [40] E. A. Dimarzio. Proper accounting of conformations of a polymer near a surface. J. Chem. Phys., 42:2101, 1965. [41] M. Doi and S. F. Edwards. The Theory of Polymer Dynamics. Oxford Press, New York, 1986. [42] A. Donev, A. L. Garcia, and B. J. Alder. Stochastic Event-Driven Molecular Dynamics. J. Comp. Phys., 227:2644, 2008.

154 [43] A. Donev, A. L. Garcia, and B. J. Alder. Stochastic Hard-Sphere Dynamics for Hydrodynamics of Non-Ideal Fluids. Phys. Rev. Lett, 101:075902, 2008. [44] A. Donev, S. Torquato, and F. H. Stillinger. Neighbor List Collision-Driven Molecular Dynamics Simulation for Nonspherical Particles: I. Algorithmic Details II. Applications to Ellipses and Ellipsoids. J. Comp. Phys., 202:737–764, 765–793, 2005. [45] K. D. Dorfman. DNA electrophoresis in microfabricated devices. Rev. Mod. Phys., 82:2903, 2010. [46] P. S. Doyle, B. Ladoux, and J. L. Viovy. Dynamics of a tethered polymer in shear flow. Phys. Rev. Lett., 84:4769, 2000. [47] B. Duenweg and A. J. C. Ladd. Lattice Boltzmann simulations of soft matter systems. Adv. Polym. Sci., 221:89, 2009. [48] E. R. Dufresne, T. M. Squires, M. P. Brenner, and D. G. Grier. Hydrodynamic coupling of two brownian spheres to a planar surface. Phys. Rev. Lett., 85:3317, 2000. [49] D. L. Ermak and J. A. McCammon. Brownian dynamics with hydrodynamic interactions. J. Chem. Phys., 69:1352, 1978. [50] P. Espanol and P. Warren. Statistical mechanics of dissipative particle dynamics. EPL, 30:191, 1995. [51] G. De Fabritiis, M. Serrano, R. Delgado-Buscalioni, and P. V. Coveney. Fluc-

155 tuating hydrodynamic modeling of fluids at the nanoscale.

Phys. Rev. E,

75:026307, 2007. [52] M. Fixman. Simulation of polymer dynamics 1. general theory. J. Chem. Phys, 69:1527, 1978. [53] M. Fixman. Construction of Langevin forces in the simulation of hydrodynamic interaction. Macromolecules, 19:1204, 1986. [54] M. Fyta, S. Melchionna, S. Succi, and E. Kaxiras. Hydrodynamic correlations in the translocation of a biopolymer through a nanopore: Theory and multiscale simulations. Phys. Rev. E, 78:036704, 2008. [55] G. Fytas, S. H. Anastasiadis, R. Seghrouchni, D. Vlassopoulos, J. Li, B. J. Factor, W. Theobald, and C. Toprakcioglu. Probing collective motions of terminally anchored polymers. Science, 274:2041, 1996. [56] M. Gauthier and G. Slater. Molecular dynamics simulation of a polymer chain translocating through a nanoscopic pore. Eur. Phys. J. E, 25:17, 2008. [57] S. Gerashchenko and V. Steinberg. Statistics of tumbling of a single polymer molecule in shear flow. Phys. Rev. Lett., 96:038304, 2006. [58] G. Giupponi, G. De Fabritiis, and P. V. Coveney. Hybrid method coupling fluctuating hydrodynamics and molecular dynamics for the simulation of macromolecules. J. Chem. Phys., 126:154903, 2007. [59] G. Gompper, T. Ihle, D. Kroll, and R. Winkler. Multi-particle collision dynam-

156 ics: A particle-based mesoscale simulation approach to the hydrodynamics of complex fluids. Adv. Polym. Sci., 221:1, 2009. [60] I. O. G¨otze, H. Noguchi, and G. Gompper. Mesoscale simulations of hydrodynamic squirmer interactions. Phys. Rev. E, 76:046705, 2007. [61] P. S. Grassia, E. J. Hinch, and L. C. Nitsche. Computer simulations of Brownian motion of complex systems. J. Fluid Mech., 282:373, 1995. [62] Y. Gratton and G. W. Slater. Molecular dynamics study of tethered polymers in shear flow. Eur. Phys. J. E, 17:455, 2005. [63] A. Groisman, M. Enzelberger, and S. R. Quake. Microfluidic memory and control devices. Science, 300:955, 2003. [64] R. D. Groot and P. B. Warren. Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation. J. Chem. Phys., 107:4423, 1997. [65] H. Hasimoto. On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past array of spheres. J. Fluid Mech., 5:317, 1959. [66] J. P. Hernandez-Ortiz, M. Chopra, S. Geier, and J. J. de Pablo. Hydrodynamic effects on the translocation rate of a polymer through a pore. J. Chem. Phys, 131:044904, 2009. [67] J. P. Hernandez-Ortiz, J. J. de Pablo, and M. D. Graham. Fast computation of many-particle hydrodynamic and electrostatic interactions in a confined geometry. Phys. Rev. Lett., 98:140602, 2007.

157 [68] J. P. Hernandez-Ortiz, H. Ma, J. J. de Pablo, and M. D. Graham. Cross-streamline migration in confined flowing polymer solutions: Theory and simulation. Phys. Fluids, 18:123101, 2006. [69] J. P. Hernandez-Ortiz, H. Ma, J. J. de Pablo, and M. D. Graham. Concentration distribution during flow of confined flowing polymer solutions at finite concentration: slit and grooved channel. Korea-Aust. Rheol. J., 20:143, 2008. [70] A. Hirsa, C. A. Lopez, M. A. Laytin, M. J. Vogel, and P. H. Steen. Lowdissipation capillary switches at small scales. Appl. Phys. Lett., 86(1):014106, 2005. [71] H. P. Hsu, K. Binder, L. I. Klushin, and A. M Skvortsov. What is the order of the two-dimensional polymer escape transition? Phys. Rev. E., 76:021108, 2007. [72] Y. Inoue, Y. Chen, and H. Ohashi. A mesoscopic simulation model for immiscible multiphase fluids. J. Comput. Phys., 201:191, 2004. [73] A. Izmitli, D. C. Schwartz, M. D. Graham, and J. J. de Pablo. The effect of hydrodynamic interactions on the dynamics of DNA translocation through pores. J. Chem. Phys, 128:085102, 2008. [74] R. M. Jendrejack, J. J. de Pablo, and M. D. Graham. Stochastic simulations of DNA in flow: Dynamics and the effects of hydrodynamic interactions. J. Chem. Phys, 116:7752, 2002. [75] R. M. Jendrejack, E. T. Dimalanta, D. C. Schwartz, M. D. Graham, and J. J. de Pablo. DNA dynamics in a microchannel. Phys. Rev. Lett., 91:038102, 2003.

158 [76] R. M. Jendrejack, M. D. Graham, and J. J. de Pablo. Hydrodynamic interactions in long chain polymers: Application of the chebyshev polynomial approximation in stochastic simulations. J. Chem. Phys, 113:2894, 2000. [77] R. M. Jendrejack, D. C. Schwartz, J. J. de Pablo, and M. D. Graham. Shearinduced migration in flowing polymer solutions: Simulation of long-chain deoxyribose nucleic acid in microchannels. J. Chem. Phys, 120:2513, 2004. [78] R. M. Jendrejack, D. C. Schwartz, M. D. Graham, and J. J. de Pablo. Effect of confinement on DNA dynamics in microfluidic devices. J. Chem. Phys, 119:1165, 2003. [79] J. L. Jones, M. Lal, J. N. Ruddock, and N. A. Spenley. Dynamics of a drop at a liquid/solid interface in simple shear fields: A mesoscopic simulation study. Faraday Discuss., 112:129, 1999. [80] R. Kapral. Multiparticle collision dynamics: Simulation of complex systems on mesoscales. Adv. Chem. Phys., 140:89, 2008. [81] R. Kekre, J. E. Butler, and A. J. C. Ladd. Comparison of lattice-Boltzmann and Brownian-dynamics simulations of polymer migration in confined flows. Phys. Rev. E, 82:011802, 2010. [82] N. Kikuchi, J. F. Ryder, C. M. Pooley, and J. M. Yeomans. Kinetics of the polymer collapse transition: The role of hydrodynamics. Phys. Rev. E, 71:061804, 2005. [83] P. R. Kramer, C. S. Peskin, and P. J. Atzberger. On the foundations of the

159 stochastic immersed boundary method. Comput. Method. Appl. M., 197:2232, 2008. [84] N. Laachi, M. Kenward, E. Yariv, and K. D. Dorfman. Force-driven transport through periodic entropy barriers. EPL, 80:50009, 2007. [85] A. J. C. Ladd and R. Verberg. Lattice-Boltzmann simulations of particle-fluid suspensions. J. Stat. Phys., 104:1191, 2001. [86] M. C. Lai and C. S. Peskin. An immersed boundary method with formal secondorder accuracy and reduced numerical viscosity. J. Comput. Phys, 160:705, 2000. [87] A. Lamura and G. Gompper. Numerical study of the flow around a cylinder using multi-particle collision dynamics. Eur. Phys. J. E, 9:477, 2002. [88] L. G. Leal and E. J. Hinch. The effect of weak Brownian rotations on particles in shear flow. J. Fluid Mech., 46:685, 1971. [89] S. H. Lee and R. Kapral. Mesoscopic description of solvent effects on polymer dynamics. J. Chem. Phys., 124:214901, 2006. [90] H. Ma and M. D. Graham. Theory of shear-induced migration in dilute polymer solutions near solid boundaries. Phys. Fluids, 17:083103, 2005. [91] A. Malevanets and R. Kapral. Mesoscopic model for solvent dynamics. J. Chem. Phys., 110:8605, 1999. [92] J. T. Mannion, C. H. Reccius, J. D. Cross, and H. G. Craighead. Conformational

160 analysis of single DNA molecules undergoing entropically induced motion in nanochannels. Biophys. J., 90:4538, 2006. [93] A. P. Markesteijn, O. B. Usta, I. Ali, A. C. Balazs, and J. M. Yeomans. Flow injection of polymers into nanopores. Soft Matter, 5:4575, 2009. [94] J. F. Marko and E. D. Siggia. Stretching DNA. Macromolecules, 28:8759, 1995. [95] A. Milchev, V. Yamakov, and K. Binder. Escape transition of a polymer chain: Phenomenological theory and monte carlo simulations. Phys. Chem. Chem. Phys., 1:2083, 1999. [96] S. T. Milner. Polymer brushes. Science, 251:905, 1991. [97] R. Mittal and G. Iaccarino. Immersed boundary methods. Annu. Rev. Fluid Mech., 37:239, 2005. [98] K. Mussawisade, M. Ripoll, R. G. Winkler, and G. Gompper. Dynamics of polymers in a particle-based mesoscopic solvent. J. Chem. Phys., 123:144905, 2005. [99] M. Muthukumar. Translocation of a confined polymer through a hole. Phys. Rev. Lett., 86:3188, 2001. [100] M. Muthukumar and A. Baumg¨artner. Effects of entropic barriers on polymer dynamics. Macromolecules, 22:1937, 1989. [101] H. D. Nguyen and C. K. Hall. Molecular Dynamics Simulations of Fibril Formation by Random Coil Peptides. Proc. Natl. Acad. Sci., USA 101:16180, 2004.

161 [102] A. E. Nkodo, J. M. Garnier, B. Tinland, H. J. Ren, C. Desruisseaux, L. C. McCormick, G. Drouin, and G. W. Slater. Diffusion coefficient of DNA molecules during free solution electrophoresis. Electrophoresis, 22:2424, 2001. [103] D. Nykypanchuk, H. H. Strey, and D. A. Hoagland. Brownian motion of DNA confined within a two-dimensional array. Science, 297:987, 2002. [104] S. B. Opps, J. M. Polson, and N. A. Risk. Discontinuous molecular dynamics simulation study of polymer collapse. J. Chem. Phys., 125:194904, 2006. ¨ [105] H. C. Ottinger. Stochastic Processes in Polymeric Fluids: Tools and Examples for Developing Simulation Algorithms. Springer, 1995. [106] P. J. Park and W. Sung. Dynamics of a polymer surmounting a potential barrier: The kramers problem for polymers. J. Chem. Phys., 111:5259, 1999. [107] S. Peng, F. Ding, B. Urbanc, S. V. Buldyrev, L. Cruz, H. E. Stanley, and N. V. Dokholyan. Discrete molecular dynamics simulations of peptide aggregation. Phys. Rev. E, 69:041908, 2004. [108] T. T. Perkins, D. E. Smith, and S. Chu. Single polymer dynamics in an elongational flow. Science, 276:2016, 1997. [109] T. T. Perkins, D. E. Smith, R. G. Larson, and S. Chu. Stretching of a single tethered polymer in a uniform-flow. Science, 268:83, 1995. [110] C. S. Peskin. Flow patterns around hear valves - numerical method. J. Comput. Phys, 10:252, 1972.

162 [111] F. Xijunand N. Phan-Thien, S. Chen, X. Wu, and T. Y. Ng. Simulating flow of DNA suspension using dissipative particle dynamics. Phys. Fluids, 18:063102, 2006. [112] I. V. Pivkin and G. E. Karniadakis. Controlling density fluctuations in wallbounded dissipative particle dynamics systems. Phys. Rev. Lett., 96:206001, 2006. [113] C. Pozrikidis. Computation of periodic Green’s functions of Stokes flow. J. Eng. M., 30:79, 1996. [114] P. Pranay, S.G. Anekal, J.P. Hernandez-Ortiz, and M.D. Graham. Pair collisions of fluid-filled elastic capsules in shear flow: Effects of membrane properties and polymer additives. Phys. Fluids, 22:123103, 2010. [115] A. Puliafitoa and K. Turitsyn. Numerical study of polymer tumbling in linear shear flows. Physica D, 211:9, 2005. [116] W. Reisner, N. B. Larsen, H. Flyvbjerg, J. O. Tegenfeldt, and A. Kristensen. Directed self-organization of single DNA molecules in a nanoslit via embedded nanopit arrays. Proc. Natl. Acad. Sci. U.S.A., 106:79, 2009. [117] M. Revenga, I. Zuniga, and P. Espanol. Boundary conditions in dissipative particle dynamics. Comput. Phys. Commun., 121-122:309, 1999. [118] M. Ripoll, M. H. Ernst, and P. Espanol. Large scale and mesoscopic hydrodynamics for dissipative particle dynamics. J. Chem. Phys., 115:7271, 2001.

163 [119] M. Ripoll, K. Mussawisade, R. G. Winkler, and G. Gompper. Low-Reynoldsnumber hydrodynamics of complex fluids by multi-particle-collision dynamics. Europhys. Lett., 68:106, 2004. [120] J. Rotne and S. Prager. Variational treatment of hydrodynamic interaction in polymers. J. Chem. Phys., 50:4831, 1969. [121] B. Roux. The calculation of the potential of mean force using computer simulations. Comput. Phys. Commun., 91:275, 1995. [122] J. N. Roux. Brownian particles at different times scales: a new derivation of the Smoluchowski equation. Phys. A, 188:526, 1992. [123] Y. Saad. Iterative Methods for Sparse Linear Systems, Second Edition. Society for Industrial and Applied Mathematics, 2 edition, April 2003. [124] D. Saintillan, E. Darve, and E. S. G. Shaqfeh. A smooth particle-mesh Ewald algorithm for Stokes suspension simulations: The sedimentation of fibers. Phys. Fluids, 17:033301, 2005. [125] R. B. Schoch, J. Han, and P. Renaud. Transport phenomena in nanofluidics. Rev. Mod. Phys., 80:839, 2008. [126] C. M. Schroeder, H. P. Babcock, E. S. G. Shaqfeh, and S. Chu. Observation of polymer conformation hysteresis in extensional flow. Science, 301:1515, 2003. [127] C. M. Schroeder, R. E. Teixeira, E. S. G. Shaqfeh, and S. Chu. Characteristic periodic motion of polymers in shear flow. Phys. Rev. Lett., 95, 2005. 018301.

164 [128] K. L. Sebastian and A. Debnath. Polymer in a double well: dynamics of translocation of short chains over a barrier. J. Phys.-Condens. Mat., 18:S283, 2006. [129] K. L. Sebastian and A. K. R. Paul. Kramers problem for a polymer in a double well. Phys. Rev. E, 62:927, 2000. [130] E. S. G. Shaqfeh. The dynamics of single-molecule DNA in flow. J. NonNewtonian Fluid Mech., 130:1, 2005. [131] N. Sharma and N. A. Patankar. Direct numerical simulation of the Brownian motion of particles by using fluctuating hydrodynamic equations. J. Comput. Phys., 201:466, 2004. [132] A. Sierou and J. F. Brady. Accelerated Stokesian Dynamics simulations. J. Fluid Mech., 448:115, 2001. [133] S. S. Skiena. The Algorithm Design Manual. Springer-Verlag, New York, 1997. [134] M. Slutsky. Diffusion in a half-space: From lord kelvin to path integrals. Am. J. of Phys., 73:308, 2005. [135] D. E. Smith, H. P. Babcock, and S. Chu. Single-polymer dynamics in steady shear flow. Science, 283:1724, 1999. [136] D. E. Smith and S. Chu. Response of flexible polymers to a sudden elongational flow. Science, 281:1335, 1998. [137] D. E. Smith, T. T. Perkins, and S. Chu. Dynamical Scaling of DNA Diffusion Coefficients. Macromolecules, 29:1372, 1996.

165 [138] S. B. Smith, L. Finzi, and C. Bustamante. Direct mechanical measurements of the elasticity of single DNA-molecules by using magnetic beads. Science, 258:1122, 1992. [139] S. W. Smith, C. K. Hall, and B. D. Freeman. Molecular Dynamics for Polymeric Fluids Using Discontinuous Potentials. J. Comp. Phys., 134:16, 1997. [140] T. M. Squires and S. R. Quake. Microfluidics: Fluid physics at the nanoliter scale. Rev. Mod. Phys., 77:977, 2005. [141] D. Stein, F. H. J. van der Heyden, W. J. A. Koopmans, and C. Dekker. Pressuredriven transport of confined DNA polymers in fluidic channels. Proc. Natl. Acad. Sci. U.S.A., 103:15853, 2006. [142] N. C. Stellwagen, C. Gelfi, and P. G. Righetti. The free solution mobility of DNA. Biopolymers, 42:687, 1997. [143] N. P. Teclemariam, V. A. Beck, E. S. G. Shaqfeh, and S. J. Muller. Dynamics of DNA polymers in post arrays: Comparison of single molecule experiments and simulations. Macromolecules, 40:3848, 2007. [144] T. Tlusty. Screening by symmetry of long-range hydrodynamic interactions of polymers confined in sheets. Macromolecules, 39:3927, 2006. [145] A. K. Tornberg and L. Greengard. A fast multipole method for the threedimensional stokes equations. J. Comput. Phys., 227:1613, 2008. [146] D. Trebotich, G. H. Miller, P. Colella, D. T. Graves, D. F. Martin, and P. O.

166 Schwartz. A Tightly Coupled Particle-Fluid Model for DNA-Laden Flows in Complex Microscale Geometries. Comp. Fluid Solid Mech., pages 1018, 2005. [147] L. N. Trefethen and D. Bau. Numerical Linear Algebra. SIAM: Society for Industrial and Applied Mathematics, 1997. [148] P. T. Underhill, J. P. Hernandez-Ortiz, and M. D. Graham. Diffusion and spatial correlations in suspensions of swimming particles. Phys. Rev. Lett., 100:248101, 2008. [149] P. T. Underhill and P. S. Doyle. On the coarse-graining of polymers into beadspring chains. J. Non-Newtonian Fluid Mech., 122:31, 2004. [150] O. B. Usta, A. J. .C Ladd, and J. E. Butler. Lattice-Boltzmann simulations of the dynamics of polymer solutions in periodic and confined geometries. J. Chem. Phys, 122:094902, 2005. [151] J. L. Viovy. Electrophoresis of DNA and other polyelectrolytes: Physical mechanisms. Rev. Mod. Phys., 72:813, 2000. [152] M. J. Vogel, P. Ehrhard, and P. H. Steen. The electroosmotic droplet switch: Countering capillarity with electrokinetics.

Proc. Natl. Acad. Sci. U.S.A.,

102:11974, 2005. [153] Y. Wang, I. Teraoka, F. Y. Hansen, . H. Peters, and O. Hassager. A Theoretical Study of the Separation Principle in Size Exclusion Chromatography. Macromolecules, 43:1651, 2010.

167 [154] N. Watari, M. Makino, N. Kikuchi, R.G. Larson, and M. Doi. Simulation of DNA motion in a microchannel using stochastic rotation dynamics. J. Chem. Phys, 126:094902, 2007. [155] S. A. Williams, J. B. Bell, and A. L. Garcia. Algorithm Refinement for Fluctuating Hydrodynamics. SIAM Multiscale Modeling and Simulation, 6:1256, 2008. [156] R. G. Winkler.

Semiflexible Polymers in Shear Flow.

Phys. Rev. Lett.,

97:128301, 2006. [157] H. Yamakawa. Transport properties of polymer chains in dilute solution: Hydrodynamic interaction. J. Chem. Phys., 53:436, 1970. [158] Y. Zhang, A. Donev, T. Weisgraber, B. J. Alder, M. D. Graham, and J. J. de Pablo. Tethered DNA dynamics in shear flow. J. Chem. Phys, 130:234902, 2009. [159] H. Zhao, A. H. G. Isfahani, L. N. Olson, and J. B. Freund. A spectral boundary integral method for flowing blood cells. J. Comput. Phys., 229:3726, 2010. [160] R. Zwanzig. Nonequilibrium Statistical Mechanics. Oxford University Press, 2001.