A hydrodynamic model of long polymer chain motion in pulsed field Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 10, № 3, 2005

A "HYDRODYNAMIC" MODEL OF LONG POLYMER CHAIN MOTION IN PULSED FIELD

V. V. Chasovskikh, L. L. Frumin Novosibirsk State University, Russia e-mail: [email protected]

S. E. PELTEK

Institute of Cytology and Genetics SD RAS, Novosibirsk, Russia

Предложена "гидродинамическая" модель движения полимерных цепей, учитывающая упругие силы натяжения цепи, призванная заменить диффузионную теорию рептаций при описании динамики длинных полимерных цепей в процессе импульсного гель-электрофореза фрагментов ДНК. Для свободно сочлененной цепи на основе статистического подхода получено термодинамическое уравнение состояния отрезка цепи в поре геля, связывающее упругую силу и плотность длины цепи в поре. Равновесие электрических сил, градиента упругих сил и сил трения для отрезка цепи в поре геля описывается нелинейным уравнением диффузионного типа. На кластере из 11 компьютеров проведены обширные численные расчеты, подвердившие способность предложенной модели "одномерной гидродинамики" описывать сложное трехмерное поведение длинных цепей ДНК в процессе импульсного гель-электрофореза. Выявлена важная роль ветвлений цепи в динамике ее движения.

Introduction

Pulsed field gel-electrophoresis (PFGE) is widely used to separate long charged DNA chains. The discovery of dependence of the DNA molecule mobility on the field period in periodic electric fields has lead to creation of such advanced separation techniques as CFGE (Crossed-Field Gel Electrophoresis) [1] and FIGE (Field Inverse Gel Electrophoresis) [2]. These techniques are based on the anomalies of velocity (mobility) that occur at certain proportions between the field period and the DNA chain size. Unlike constant field electrophoresis that permits to separate only rather small molecules (up to 30... 40 Kilo base pairs, Kbp), periodic or nonuniform electric fields can separate chains of many hundreds and thousands Kbp.

At the same time, despite a fifteen-year history of research and its special bearing on biological applications, the physical causes of such mobility anomalies, as "antiresonance" (a deep minimum of the mobility) and band inversion, have not been comprehensively studied so far. The detailed experiments by Sabanayagam and Holzwarth [3] who measured instant velocities of the DNA fractions during positive and negative field pulses practically have not received interpretation. A recent research of nonlinear electrophoresis [4] and discovery of the non-monotonous dependence of the nonlinear mobility on molecular sizes have increased a

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.

number of unsolved problems and have made the modeling of pulsed field gel electrophoresis especially urgent.

Numerous theoretical models are described in the literature that were used for the description of the DNA motion in gel, both in constant and periodic electrical fields. The volume restriction of this article prevents us from the enumeration of all published works and existing approaches. Therefore, we limit ourselves to their outline. The rather complete recent review of existing models and theories can be found in the recent overview by Viovy [5].

All these models can be divided into the three groups:

— direct imitation modeling of the polymer chain motion in medium with obstacles. The simulations by Deutsch [6] belong in this class;

— diffusion models or models of reptation in gel "tubes". These are the biased reptation model (BRM) [7, 8] and the more advanced biased reptation model with fluctuation (BRF) [9, 10], and also models of the repton diffusion by Duke etc. [11-13];

— the models that take into account elastic "entropic" forces of the polymer chain stretching. There are only few of these models. One can relate hither the "Lakes and Straits" model by Zimm [14, 15] and the works by Slater and co-authors [16].

We will not consider the direct imitation modeling, because it requires a large number of poorly known parameters describing elastic properties of the polymer chain and its interaction with gel and Brownian motion. Complexity of such approach and the essential arbitrariness of parameter settings significantly limit its usefulness even for the motion of rather short chains.

Diffusion models, similar to the BRF and lattice models, employ far more terse chain description. They are based on the reptation mechanism by de Gennes [17] and the "tube" conception by Doi and Edwards [18]. The motion of a long polymer chain in these models is reduced to the diffusion- and migration-like motion in gas of conformation "defects".

Diffusion models have made it possible to give quantitative interpretation to the motion of small DNA fragments in constant electric fields. These models are proved in numerous electrophoretic experiments and give an explanation to many steady-state electrophoresis data and such effects as band inversion in weak constant fields, for example. The data on the effective DNA charge density and characteristic gel pore sizes, obtained in the scope of these models, is important in itself.

Unfortunately, the diffusion models can not describe a deep (almost tenfold) minimum of mobility in a session of FIGE. As concerns long chains and periodic electric fields, linear diffusion models are not applicable here, in our opinion, as they do not take into account elastic forces of a chain tension.

The difficulties of diffusion models become obvious when we consider a U-shaped conformation of the charged polymer chain in an electric field. The Boltzmann distribution of reptons or defects leads to an (exponential) energy barrier for U-shaped conformations. This barrier prohibits the motion of defects and of the whole chain. This difficulty is noticed in [5], but as a problem of the Monte-Carlo method. There were attempts to overcome the diffusion model drawbacks [13, 14] with the help of the approach dubbed the "non-local Monte-Carlo" method. This approach allows for repton "jumps" of any distances, but we believe it not sufficiently proved.

The validity criterion of diffusion models is small change in the density of defects (repton) found on a chain of a length L placed in an electric field E. Adopting the Boltzmann distribution for the density of defects, this criterion can be written down as

qEL/T < 1,

(1)

where q is the charge of a defect and T is the temperature of medium in energy units. The estimates of criterion (1) show that diffusion models for typical electrophoretic fields can describe only rather short DNA fragments with a size of less than several hundreds of Kuhn lengths.

Thus, a major element of the PFGE description of long DNA chains is the contribution of elastic tension forces, which are frequently called "entropic forces" after Doi and Edwards [18]. Slater and co-authors [16] were possibly one of the first who suggested taking into account the tension forces of a polymer chain, with reference to DNA molecules. In the above-mentioned "Lakes and Straits Model" by Zimm [14, 15], these forces are also included but semi-empirically, because the Gaussian chain model is not working in the case of a short polymer chain segment in a gel pore. In this connection, we believe that the freely jointed chain model is more adequate. It is certainly the case for not very high percentage of agarose gel, when one pore makes room for a few Kuhn lengths of the chain.

In this work, we consistently use the freely jointed chain model. To describe a chain segment, we calculate the statistical integral and derive the equation connecting the elastic force and the density of length in a gel pore. This equation can be treated as a thermodynamic "equation of state" for a chain segment in the pore. Use of the equation of state is characteristic for the hydrodynamic approach. Therefore, one can speak about a hydrodynamic model of the DNA motion. It is known that the diffusion approximation is quite helpful for the description of the motion of matter under small pressure forces. However, with the pressure forces growing, the hydrodynamic approach becomes necessary, which includes the equation of state connecting the pressure forces and the matter density. By analogy, the motion of a short polymer chain under the action of small elastic forces can be described with the diffusion approximation. But the essentially different approach similar to the hydrodynamic approach of fluid dynamics is necessary when the elastic forces become great.

The purpose of this work is to construct a model correctly taking into account the elastic forces and permitting to describe key experiments on gel electrophoresis of long DNA molecules.

1. The statistical mechanics of a freely jointed polymer chain segment in a gel pore

1.1. Freely jointed polymer chain

In the rest of the article, we will consider a freely jointed chain consisting of rigid links with a characteristic size of b. With reference to a DNA molecule the size of a link is adopted to equal the Kuhn length. We suppose that the chain length, L, considerably exceeds the link size, b. The polymer chain is placed in gel, which is a set of pores (cells). The gel pores will be characterized by their root-mean-square dimension h > b and by the ratio

The polymer DNA chain is uniformly charged with a linear density of charge, x The charge of one segment is given by the expression

For the 1-percent agarose gel, the parameters h and b can be adopted to equal approximately ~ 300 ^m and 100 ^m, respectively. This gives a value of about 3 for the m parameter. As

m = h/b.

(2)

q = xb-

(3)

an estimation of the x magnitude, we will adopt a value of 250... 330 elementary charges per 1 ^m. Let us estimate ratio e of the electrical energy per segment, qEb, to its heat energy. In fields of 1... 10 V/cm this ratio will make e = qEb/T ~ 10-2 ... 10-1. We see that for all the field strengths, usually used for electrophoresis, e < 1. The small value of e enables us to develop perturbation techniques with this parameter.

For the description of a chain state, we use a ratio w that gives the density of length for a chain in an individual gel pore:

w = l/h, (4)

where I is the length of a chain segment trapped in a given gel pore. This length is just the very matter, which flows from one pore to another and leads to the overall chain motion. For a polymer chain, the density of length plays the same role as the density of mass plays in fluid dynamic problems. It is convenient to define a number n of segments trapped in one pore

n = l/b. (5)

Using ratios (2), (4), and (5), we present the density of length as

w = n/m. (6)

The ultimate value for the density of length w = 1 corresponds to the completely overstretched chain (n = m).

In the absence of an electric field, the DNA chain has an equilibrium conformation of a stochastic coil with the characteristic size R2 = bL. For the following, it is of interest to determine the equilibrium density of length, w0, in a gel pore in the absence of a field. For the length I of a freely-jointed chain segment, the root-mean-square size of which is equal to the pore size, one can write down

h2 = bl. (7)

This expression, which can be presented as

m2 = n, (8)

is exact for the free-jointed chain [19]. Its comparison with the expression for a more realistic (for the DNA case) model of an elastic chain by Blesler — Frenkel [20] shows that the validity range of the freely-jointed chain model over the size of one gel pore is given by the inequality:

n > 1. (9)

This inequality limits agarose gel concentration to a value of about < 1 %.

Substituting expression (8) into equation (6), we obtain, for the equilibrium density of length in the absence of a field:

w0 = m. (10)

In the absence of a field, chain has an equilibrium conformation of a stochastic coil and no elastic stress arises in the chain. After applying an electric field, the chain begins to flow from one pore to another, the non-uniform distribution of density of length arises, and elastic stress does appear. The reason of the elastic stress is Brownian motion of the chain. In the case of freely-jointed chain this is Brownian motion of the chain links.

The dynamic of change of a chain configuration is slow in comparison with Brownian fluctuations of its links. Therefore we fix for all cases that each chain segment found in a

gel pore has reached its local thermodynamic equilibrium (LTE) with its environment. The LTE assumption appears to be reliably established if one takes into account the characteristic times of Brownian motion for DNA segments. It permits to apply the statistical mechanics methods to examine the chain state and to calculate statistical integral for a chain segment occupying one gel pore.

Unfortunately, the statistical integral calculation cannot take into account both elastic and electric forces because of great mathematical difficulties. One can carry out the statistical integral calculation for the case when there is a noticeable elastic force and the electric force acting on the chain segment in the pore is small compared to the elastic force. This is the case within the internal range of the chain, where the tension force usually considerably exceeds the electric force. As concerns endpoints of the chain and its arms, one can neglect the elastic force, and the electric force determines the statistics of these parts of the chain.

1.2. Statistical integral of a freely jointed polymer chain segment in the elastic force case

Consider the statistical integral Z of a chain segment stretched by a force f with no electric field contribution to the segment energy. The chain segment is represented by n vectors bk and corresponds to one gel pore. We choose the z-axis of the Cartesian coordinate system along this gel pore. We denote by 9k the angle between vector bk and the z-axis, and by 0k the azimuth angle with respect to the x-axis of the Cartesian coordinate system (x,y,z).

For the statistical integral, we write down the expression

Z = y"exp^- fbk/T^Jdtt, (11)

where T is the temperature in energy units, dQ = nk[sin(0k)d9kd<fik/4n] , and fbk = fbcos(0k).

The index k runs from 1 to n in the sum and in the nk product.

k

The statistical integral (11) permits to find the average of (cos(0)), which determines the relative length w of a chain segment in the pore. For (cos(0)), expression (11) gives

(cos,)) = M ■ (12)

where the dimensionless force ft is introduced:

ft = fb/T. (13)

The direct calculation of the statistical integral gives the expression:

Z = [sinh(ft )/ft]n. (14)

Using equation (12) for (cos(0)), we obtain

(cos(0)) = coth(ft) - 1/ft = A(ft). (15)

Here, A(ft) is the well-known Langevin function. This result is also well known and leads to the classical equation connecting the elastic force and the chain extension A/ = nb(cos(0)):

A/ = nb A(ß ).

(16)

This very equation is commonly used to describe the elastic "entropic" forces in a freely jointed chain consisting of the fixed number n of rigid links. We consider the equation (16) to be inadequate for the description of the chain segment in the gel pore, because the number n of links in the pore is not constant and depends on applied forces. Indeed, equation (16) belongs in the statistics of the system with the fixed number of particles (the chain with the fixed number of links), while we are interested in the system with the variable number of particles and with the volume fixed: a chain segment in one gel pore. It is necessary to find out the dependence of the density of length, w, for a chain segment in the presence of the elastic force f in one gel pore. We find the density of length from the following equation, which determines the root-mean-square size h of the pore

h2

+

(17)

Here, the sizes xk, yk, zk are the Cartesian projections of a k-th link bk of the chain in the pore, and the angular brackets (•) mean averaging over all statistically possible configurations of links in the chain. In the spherical system of coordinates we have relationships for

the coordinates of the link

Xk = b sin(0fc) sin(0fc), yk = b sin(0fc) cos(0fc), Zk = b cos(0fc). Now, equation (17) can be presented in the form

(18)

h2 = b2

^T cos(6>k)

+

J^sin^k )cos(0k )

+

J^sin^k) sin(0k)

. (19)

Taking into account that average values (cos(0k)) = (sin(0k)) = 0, expression (19) can be reduced to:

h2/b2

m

n +

^cos(6>k )

E<cos2(^k)> .

Then, after simple transformations we obtain

m2 = n + ^ (cos(9i))(cos(9k ))•

(20)

i#k

At last, taking into account that (cos(0k)) is given by expression (15) and does not depend on its index k, equation (20) can be transformed to the more convenient form for calculation:

m

n +(n2 - n)A2(ß).

(21)

Equation (21) implicitly determines the density of length of the chain segment, w = n/m, as the function of force f. We transform equation (21) to the form

A(ß) = [(m2 - n)/(n2 - n)]1/2.

(22)

Using equation (6) for the density of length w and equation (10) for the equilibrium density of length wo in the equilibrium state of a stochastic coil, equation (20) can be written down as

A(ß) = [(wo/w - 1)/(wwo - 1)]1/2.

(23)

2

2

2

2

We have obtained, in the implicit form, the dependence of the elastic force on the density. This dependence can be called the thermodynamic equation of state for the chain by analogy with the dependence of the pressure on the density for gas. With the inverse Langevin function V the equation of state can be presented in the form

In fig. 1, the dependence of the elastic force on the density of length in the pore is plotted. Consider in more detail the asymptotic of the equation of state for the cases of strongly and weakly stretched chain segments. For a weak tension, using the asymptotic of the Langevin function A(P) ^ P/3 with P ^ 0, we obtain the equation

As we expect, the tension force tends to zero when the density tends to its equilibrium value w0 in the absence of a field.

In the case of the maximum stretching of the chain, w ^1, there occurs the indefinitely great force of chain tension. Uses the asymptotic A(P) ^ 1 — 1/P with P ^ to, we obtain

The latter equation is quite opposite to the equation of state of ideal gas in that the elastic forces do not decrease with density w decreasing. On the contrary, these forces increase. We use the equation of state, (24), to calculate the elastic forces in the internal overstretched chain links, when describing the chain motion in section 3. The equation of state permits to determine the chain tension force for the given density of length w in the pore. Notice also that the entropic elastic force f represents an analog of one-dimensional "pressure" force, which is always parallel to the generatrix of the gel pore, in which a polymer chain segment resides.

Taking into account that the antiresonance is more pronounced for weaker fields and decreases when the field increases, the explanation for the deep minimum of the mobility in

f = (T/b)V{[(wo/w - 1)/(wwo - 1)]1/2}.

(24)

f = 3(T/b)[(wo/w - 1)/(wwo - 1)]1/2.

(25)

f « (T/b)/(w - 1).

(26)

/[KT13 N]

10

5

Fig. 1. Equation of state: dependence of the elastic force f on the density of length w in the pore, w0 = 3.

antiresonance does not mandate including the effect of an electric field into the equation of state for the chain. Notice that the elastic force at the chain endpoints is small and the electric field must be taken into account. The following section of this paper deals with this matter.

1.3. Statistics of the ends of charged freely jointed polymer chain in electric field

Now, we consider free end segments of the chain, where the elastic forces tend to zero. The electric field, acting on the charged segments of chain endpoints, aligns these segments along the field. Let this (endpoint) chain segment be represented by n vectors bk. We direct the z-axis of a Cartesian coordinate system along the electric field, E. The value 0k is the angle between vector bk and the z-axis.

The potential energy of the first link in the electric field with respect to this link's beginning is equal to U = —Eqbcos($i)/2. For the second link we arrive at U2 = —Eqb(cos(0i) + cos(02)/2). The potential energy of a k link is given by the sum

Uk = —Eqb

J^cos(^) + cos(0k )/2

,i<k

(27)

The complete potential energy of the chain for the given array of link angles © =

02,..., 0n} is equal to the sum U = Uk

k

U = Eqb[(n — 1/2) cos(0i) + ... + (n — k + 1/2) cos(0k) + ... + (1/2) cos(0n)]. (28)

Omitting the constant factor 1/2n, one can present the statistical integral Z for the chain segment in the form

Z = J exp(—U/T)dQ, (29)

where dQ = k { sin(0k)d0k}.

Integral (29) can be easily found. For the following, however, instead of the statistical integral, it is more convenient to calculate the many-parametric statistical integral, S, which we define with the help of an n-dimensional vector e = {ek}:

•(e) = exp

Y,ek(n — k + 1/2) cos(0k)

dQ. (30)

The calculation of many-parametric statistical integral (30) leads to the expression:

S(e) = J] { sin h[ek(n — k + 1/2)]/[ek(n — k + 1/2)]}. (31)

k

All statistical properties of the chain segment are described by this many-parametric statistical integral. For example, the statistical integral itself can be derived from many-parametric statistical integral S(e), if we assign ek = e = Eqb/T.

The many-parametric statistical integral permits to find the averaged (cos(0k)):

J_d ln[S(e)]

n — k + 1/2 dek

(cos(0k)) =_,.,, /0 r =Eqb/T. (32)

From the equations (31) and (32) it follows that

(cos(tffc)) = A[e(n - k + 1/2)], (33)

where A is the Langevin function. Expression (33) shows that the averaged cosine depends on an ordinal number of a link. The "higher" with respect to the field the link is found, the more orderly it is oriented along the field. Thus the bottom link (k = n, (cos(0n)) = A(e/2)) experiences the weakest orienting exertion by the field.

Let us proceed now with the calculation of the density of length at chain endpoints in the presence of the field. We will use equations (17)-(20) and transform the equation (20) to the form more convenient for calculation:

m2 = n

+ )»2 - ^«cos(0fc))\ . (34)

V k k J

From equation (33), we obtain for <cos(0k))

m2 = n + j^A[e(n - k + 1/2)] | -J] A2[e(n - k + 1/2)]. (35)

k

This equation implicitly determines the density of length w = n/m at the chain endpoints as a function of an electric field strength E. As one sees from (35), in the strong field limit, using the asymptotic of the Langevin function, A[e] ^ 1 at e ^ we obtain the completely overstretched chain segment with number of links n = m.

In the weak field limit we have n = m2. It means that in the weak fields (e ^ 1) one can adopt as a boundary condition for density of length:

wr = w0. (36)

2. Equilibrium of forces and discrete model of polymer chain motion

2.1. "One-dimensional Hydrodynanics" of the polymer chain motion

Consider the equilibrium of forces applied to a polymer chain segment. We will take into account that the elastic force f plays a role of the one-dimensional pressure force. As well as in the case of hydrostatics, where specific forces applied to a unit of mass counterbalance the gradient of pressure, the one-dimensional equations of the equilibrium of forces applied to the chain segment can be written down in the form

f = -K = -(Ke + KF ). (37)

Here f is the elastic entropic force, s is the coordinate along a gel tube, in which the chain segment is placed, K is the component of external forces in the pore direction, referred to a unit of pore length. The force K is a sum of the specific electric force component, KE, and the specific friction force, KF, acting on the chain from the medium.

To calculate the specific electric force, we take a gel pore with a size of h = mb, which contains n Kuhn segments of the polymer chain. The electric force component in the pore direction, divided by the pore size, is equal to

KE = nbxE cos(9)/h = (n/m)xE = wxE cos(9), (38)

where 9 is the angle between the pore direction and the electric field direction.

The specific friction force, KF, applied to the chain segment from the gel, is not the function of the density of length, but is proportional to its flux J. It is determined by the expression

Kf = -J, (39)

where £ is the effective friction coefficient per unit of the pore length, which takes into account both the friction against gel and the Stokes friction.

Combining expressions (37)-(39), we obtain the equation for the flux of length in the pore.

df

J = pwxE cos(9) + p^-. (40)

ds

Here, factor p = 1/£ is similar to the mobility. The existence of the flux of length, J, leads to a change in the length density and results in the overall motion of the polymer chain.

Finally, we consider the chain length conservation law in the form of one-dimensional equation

dw dJ . A .

u + T. = ° (41)

The polymer chain motion in gel is completely determined by the joint solution of equations

(40) and (41). Substituting (40) in equation (41), we obtain the equation

dw d [w cos(9)] d 2f

li + px^^s^ + PTS-2 =°. <42)

To describe the chain motion one should use equation (42) or the system of equations (40),

(41), together with the equation of state for the chain, (24).

Equation (42) resembles the one-dimensional equation of defect diffusion by de Gennes [19]. The difference is that (42) takes into account the elastic forces in the polymer chain and represents hydrodynamic description of the chain motion in gel, although formally belongs in the class of one-dimensional (nonlinear) diffusion equations.

The likeness of equation (42) and de Gennes' diffusion equation of gas of defects explains the success some attempts had in expanding the scope of diffusion models beyond the limits of their validity, for example, to long polymer chains and large electric fields.

It is important to emphasize the nonlinear character of equation (42). As we have noted, recent research of the nonlinear mobility of DNA molecules [4] has not given the interpretation within the framework of linear diffusion models. Nonlinear equation (42) leaves an opportunity for the explanation of these effects.

2.2. Discrete model of polymer chain motion

The differential equations obtained in section 1.3 have the formal character so far, since the elastic entropic force is determined only for sufficiently bulky chain segments of a few Kuhn

lengths. Therefore the discrete analog of equations (40) and (41) appears to be richer in content from the physical point of view. The natural quantization size is the root-mean-square size h of a gel pore.

Let the polymer chain occupy N pores at some moment. We denote the density of length in a k-th pore by wk and write down the discrete analog of mass transfer equation (41) in the form of an explicit scheme

(wk+1 - wk)/t = - (Jk+1/2 - Jk-1/2) /h. (43)

Here, top index i is used to enumerate time steps (time layers). Fluxes J refer to the pore middles and are enumerated with half-integer indices. The time step interval is denoted by t. Solving equation (43) for the unknown value of the density of length at an (i + 1) time step, we arrive at the first equation of our discrete model

wk+1 = wk - h Jk+i/2 - Jk-1/2) • (44)

Each pore is characterized by vector r, the module of which is equal to the root-mean-square size h of the pore. The direction of vector r determines the spatial orientation of the pore. The sum of vectors r gives an end-to-end chain vector.

Lowering the top index, i, we write down the difference approximation of equation (40) for the fluxes in the form

Jk+1 /2 = cos 9k+1/2(wk+1 + wk)/2 - Mfk+1 - fk)/h,

Jk-1/2 = cos 9k-1/2(wk + wk-1)/2 - Mfk - fk-1)/h- (45)

Here 9 is the angle between vector r and an electric field E.

The boundary conditions for the system of equations (44), (45) are given by equation (36).

Now, let us describe the chain motion calculation procedure. In the early moment, before a field is turned on, the chain represents a stochastic coil that occupies N intercommunicating pores randomly oriented in space. All the pores occupied by the chain are enumerated from the "tail" to the "head" of the polymer chain.

At each time step, the system of equations (44) and (45) is solved with account taken for the boundary condition (36). For the consecutive pores occupied by the chain, the values of density and flux are calculated. To solve the equation of state (24), the inverse Langevin function was tabulated for 200 values of the density of length in a range from 1 to wo. In between the selected values, the linear approximation was used.

During the solving of the system of equations (44), (45), we constantly watched the fluxes of length into both the first and the N-th endpoint pores. The density of length in the endpoint pores, wr, is given by the boundary condition (36). The flux of length into the endpoint pore leads to the generation of an "excessive" length there. If this excessive length, which flowed in into the endpoint gel pore, exceeded the boundary value hwr for the length of the endpoint chain segment, a new occupied pore was added with the density equal to wr. Excessive length can be negative as well, thus leading to a decrease in the number of pores occupied by chain. Such algorithm keeps constant the total (contour) chain length L, although the number N of pores that the chain occupies at various moments during calculation can vary.

When the chain endpoints enter into and exit from the gel pores in their run, the chain moves in space. In the course of modeling, the special attention was paid to the procedure of selecting vector r for a new pore. When the field is weak or there is no field at all, this vector

should be generated with an isotropic orientation in space. The strong electric field stretches chain segments, and the direction of the new pore occupied by the chain becomes anisotropic in space. For the probability density W of the angular distribution, ref. [9] suggests the use of the angular distribution for a rigid rod of size le with the probability density

The statistics of configurations for endpoint segments of the freely jointed chain differs from the statistics for a rigid rod. But due to its simplicity, the model of the rigid rod with a fitted effective length is quite suitable for choosing an angle for the new pore.

The explicit scheme of (44), (45) has the second order of accuracy by spatial variable and the first order by time variable. This scheme is conservative, i.e. keeps the length of the chain constant during the calculation, but is not absolutely steady. The stability of this chain depends not only on a proportion of time and spatial steps, but also on a value of the density of length, since the equations are nonlinear. The criterion of stability of the explicit scheme for the diffusion equation has the form [21]:

where D is the diffusion coefficient. Due to non-linearity of our problem, its effective diffusion df

coefficient D = — — is a non-trivial function of the density of length. During calculation dw

the criterion (47) was constantly reviewed and the scheme parameters (time step t and spatial step h) were duly adjusted for maintaining scheme stability. It is interesting to note, that while the "common" hydrodynamics requires a denser grid for larger pressure gradients, in our case the stability of calculations improves when taking a sparser grid. This fact is explained by the fundamental difference of the equation of state for a polymer chain segment (24) from the equation of state for ideal gas.

Despite the poor stability of the explicit scheme, the use of steadier implicit schemes has proved to be impossible for the following reason. Tube leakages, or "hernias", which accompany the chain motion, require taking into account the interaction of polymer chain segments trapped in one pore. This non-local interaction of chain segments with one another prohibits the use of the implicit schemes. The chain hernias are examined in the following section.

2.3. Tube leakages or, "Hernias"

The discrete equations of motion and boundary conditions described above permit to simulate the polymer chain motion using only one empirical parameter, either friction coefficient £ or its inverse parameter analogous to the mobility. At the same time, the proposed one-dimensional hydrodynamic description of the chain motion is not complete, because it does not take into account tube leakages, or "hernias".

The tube leakages in our model arise naturally when a significant overflow of the density of length, w > wr, occurs in some pore occupied with a chain segment. With this pore overflow, no elastic forces arise in the freely jointed chain model, but the part of the chain begins to protrude (herniate) into the adjoining pore. This protrusion leads to the formation of U-turn on the chain, which we call "hernia". The calculation scheme admits the formation of multiple hernias, when a new hernia is formed from the existing one.

Consider in more detail the modeling of hernia growth. When the condition of hernia forming (w > 2wr) in an i-th vertex of the chain is met, one new pore and two new vertices are added.

W - exp[xEl2 cos(0)/T].

(46)

t < h2/2D,

(47)

One new vertex is formed at the end of the new pore (this is the hernia vertex) and another one in the beginning of the pore. In the hernia vertex, the density is adopted to equal the boundary density wr.

The mechanism of hernia destruction (dissolving) is realized in a similar fashion. In its motion, the hernia can occupy multiple pores; secondary hernias and even hernia bunches can be formed. The interaction between chain segments occupying a common pore becomes important. Really, simple expression (39) connecting the specific friction force and the length flux does not work now, because the hernia occupies the pore and there are at least two polymer chain segments, which interact with each other.

Consider the two chain segments occupying a common pore. The length fluxes in each of these segments are denoted by Ji and J2. The additional friction force Ki2, acting on the first segment from the second, should depend on the relative velocity of the two segments. Because these segments are oriented opposite to each other (r1 = — r2), this force is proportional to the sum of the length fluxes

K12 = — £c (Ji + J2), (48)

where is the effective coefficient of friction between chain segments per unit of the pore length. The equation of the force equilibrium for these segments can be written down in the form

Ji + £c (Ji + J2) = Kei + f,

J2 + £c (Ji + J2) = ke2 + f. (49)

Here, KEi and KE2 are the specific electric forces acting on the respective chain segments from the field. Let the right parts of equations (49) be Ki and K2, respectively. With these designations, the joint solution of the equations for the Ji and J2 fluxes is given by the expressions

Ji = ^Ki — ^c (Ki + K2)/(£ + £c ),

J2 = ^K2 — ^C (Ki + K2)/(£ + £c ). (50)

In practical calculation, we adopted that friction coefficients for the chain against a gel fiber and for the chain against itself are not very different: « We used the discrete analogue of expressions (50) for calculation of the chain motion in hernias.

Expressions (50) represent non-local interaction, because the segments marked by the indices 1 and 2, respectively, are not the adjacent chain segments. They interact with each other because they happened to hit the same pore during the hernia formation.

The non-local interaction in hernias complicates "sliding" of segments against each other and leads to slowing down of "dissolving" of hernias. This effect slows down the chain motion and, in the end, results in stabilization of U-shaped conformations, which are formed during change of an electrical field polarity. The extensive numerical experiments have shown, that the internal friction of a chain in pores with hernias is responsible for the deep minimum of mobility at "antiresonance".

The noticeable role in the hernia formation is played by an effect of the excluded volume [21]. Numerical modeling revealed that an alternating electric field forms up to 7 hernias in some pores simultaneously. Taking into account the excluded volume effect becomes necessary. This effect consists in decreasing volume available for new chain segments in a pore occupied by the

other segments. This decreases the effective boundary density wreff. The effect is accompanied by decreasing threshold of the secondary hernia formation. The secondary hernias are formed earlier and their total amount increases. We simulated the effect of excluded volume as a decrease in the effective boundary density, with the help of an approximation

wreff =1 + (wr — 1)/ka, (51)

where k is the number of segments occupying the given pore. A critical parameter, a, was assumed to equal 0.5. The calculation showed that the excluded volume effect rather weakly affects the minimum velocity in antiresonance, but leads to noticeable (up to 20 %) shift of the antiresonance period.

3. The calculation results

We carried out the calculation on the basis of the above hydrodynamic model, which permitted to quantitatively interpret some of the most interesting experiments on PFGE of DNA long chains in agarose gel. The calculation model served as a microscope-like tool and permitted to observe the behavior of a charged long polymer chain in periodic fields.

The principal result obtained is that the calculation confirmed the possibility of deep minimums of the mobility in moderate (2... 6 V/cm) electric fields. The calculations were performed with the "classical" realization of FIGE, when the duration T+ of a direct (positive) pulse of the field is 3 times the duration T- of an inverse (negative) pulse (Rt = T+/T- = 3). The amplitudes of positive E+ and negative E_ field pulses were identical (RE = E+/E_ = 1).

During calculation, we observed the motion of an individual DNA chain over many field periods. Under assumption of process ergodicity, the time-averaged chain motion must coincide with the ensemble-averaged motion for a large ensemble of chains, which is observed in real experiments.

The qualitative description of the chain motion during numerical experiments, on the first sight, is not much different with the scheme suggested long ago in [18], which assigned the special role to a U-shaped chain conformations. The essential difference consists in the extremely important role of hernias for the chain motion.

The chain practically always travels over a path beaten by one of hernias. During the chain motion, hernias always arise and disappear. Those oriented in the field direction grow quickly, while those oriented at angle gradually "dissolve". Such competition results in the strong orientation of a chain in the field direction. The numerical experiments have shown, that even if the occupation of new pores were isotropic, after this competition of hernias, the chain becomes noticeably oriented in the field.

In constant fields, the motion of a long chain is rather uniform, but accompanied with noticeable fluctuations, and quite often U-shaped conformations are formed. The "head" and "tail" of a chain frequently change places. It was observed often in the real experiment. In periodic fields, the U-shaped conformations practically always arise when the field direction flips. These conformations are responsible for "antiresonance". However, a deep minimum of the mobility in antiresonance can be obtained only if we take into account the internal friction of a herniated chain. Without this friction, hernias quickly disappear, and we observe the high velocity of the chain motion and a shallow minimum in "antiresonance".

The phases of the chain motion are shown in fig. 2.

{ÄTC^O,^^ ■Sj i=0s

i= 1.3 5 s

i=2.7s

i=4.05s

i=5.4s

i=6.75s

i=8.1s

J_ -— - i=9.45s

J— '-"—'—^ i=10.8s

10

20

30

E, V/cm

-1--r>

-4 +4

40

X, mkM

Fig. 2. The phases of the DNA chain motion in a periodic electric field. The calculation corresponds to a chain size of 340 Kbp. The field strength is E+ = 4 V/cm, Re = E+/E- = 1, Rt = T+/T- = 3. Duration of one cycle is 10.8 s. The hernias of the chain are shown in thick lines.

In "antiresonance" the chain oscillates in a fixed place for the most part of periods. At the end of each positive pulse of a field it acquires the distinct U-shaped configuration. The greater the time a chain spends in this U-shaped trap, the deeper the minimum of the mobility at "antiresonance". We emphasize that practically at all the phases the chain motion is determined by the growth of and the competition between hernias, and the non-local interaction in the latter results in stabilization of U-shaped conformations.

In fig. 3, the calculation of the DNA velocity as a function of the positive pulse period of a periodic field is shown for the standard FIGE case. The calculation is performed for a DNA size of 340 Kbp. The velocity was calculated as a ratio x/T of the chain center-of-mass displacement, x, to a time interval of T.

The triangles in the figure denote the center-of-mass velocity of a molecule during a positive (direct) field pulse, the squares correspond to the velocity during a negative (inverse) pulse, and circles mark the resulting curve of the velocity averaged over the field period.

Fig. 3 demonstrates a distinct, deep minimum of the chain's center-of-mass velocity. The velocity behavior during both positive and negative field pulses corresponds well to unique experiments performed in ref. [3]. The values of characteristic drift velocities and the time corresponding to the minimum mobility also coincide rather well with the data of these experiments. The results shown in the figure were obtained by averaging the velocity of one molecule over the long time intervals, from 1000 to 5000 s. The averaging time was longer for points near the minimum and shorter for the long periods of a field. The figure shows that for longer periods the average velocity is almost constant and does not depend on the field period.

In fig. 4, the drift velocity dependence on the period of a positive cycle for electric field strengths of 2, 4, and 6 V/cm is presented.

X/t, mkM/s

5 4 3 2 1

2 4 6 8 101 2 4 6

T+, s

Fig. 3. The velocity of a 340 Kbp DNA fragment as a function of pulse time T+. E+ = 4 V/cm, E+ = E- (Re = 1), Rt = T+ /T- = 3.

X!t, mkM/s

Fig. 4. The velocity of a 340 Kbp DNA fragment as a function of positive pulse time T+ and electric field strength E. RE = 1, RT = 3.

M

a

- 0

-5

- -10

- -15

"T1

-5

-10

-15

10

20

10

_L

20 t, S

Fig. 5. Dynamics of change of the DNA chain instant velocity during one cycle of a periodic field. E+ = 4 V/cm, E- = E+ (Re = 1), Rt = 3. The dashed lines below and above the abscissa axis represent average velocities for the positive and negative periods of a field pulse, respectively. The dashed line, next to zero, describes the total velocity averaged over the field period.

The figure demonstrates a deep minimum of the velocity in weak and moderate fields (2... 6 V/cm). Comparison with the ref. [3] shows the exclusive qualities of the suggested model, which permits to qualitatively and quantitatively describe the antiresonance phenomenon.

In ref. [3], also the instant drift velocity was measured during the field period. It develops a characteristic peak just after change of field polarity and the following wide minimum. It is of interest to compare these results with ours. In fig. 5, the corresponding instant velocities averaged over the large number of the field periods are presented for a 340 Kbp DNA chain.

The sharp center-of-mass velocity peaks that come immediately after changing a field polarity are caused by the fast dynamics of the chain motion, which is finished with the formation of new hernias in the field direction and a decrease in the velocity. The peak of the instant velocity immediately after changing a field polarity is caused by the drastic shrinkage of overstretched "arms" and the redistribution of the density of length along a chain. The elastic entropic force plays a definitive role in this dynamics.

In the Figure, the three plots a, b and c are shown, which correspond to the field periods shorter than resonant (a), next to resonant (b), and longer than resonant (c), respectively. For the field periods near to the resonance time, the highest instant velocity peak occurs during a negative field pulse and the wide instant velocity minimum occurs during a positive field pulse (fig. 5, c). Comparison with the experiments that were performed in ref. [3] showed that our calculations agree with the experimental data in this case as well.

4. Discussion

The presented model permitted to quantitatively describe the deep velocity minimum observed in standard FIGE sessions, velocity dependence on electric field amplitude, as well as the fast dynamics of instant velocities after changing a field polarity. The overall calculation results correspond well to experimental data and prove the validity of the suggested model.

Some restrictions of the suggested model are connected with a large volume of calculations and, consequently, with large calculation time. The presented results were obtained with the cluster of 11 360 MHz IBM compatible computers. The calculations were carried out for several months. Because of the large time the calculations for super long molecules (longer than 1 Mbp) were not complete and are not presented in this paper. As concerns shorter chains, the quantization errors become essential here.

One must understand that the accuracy of the suggested approach that use the freely jointed chain model is limited by big pores and low percentage of agarose content in gel. The model gives more accurate results for weaker fields. In strong fields, the calculation with considerably shorter time steps is required.

Conclusion

The basic results of work are as follows.

The one-dimensional "hydrodynamic" model of the long polymer chain motion that takes into account the tension forces in a charged chain is suggested. On the basis of the statistical approach, the statistical integral is calculated for the freely jointed chain model and the thermodynamic equation of state for the chain segment in the pore is obtained that connects the elastic force of chain tension and the density of length of the chains in the gel pore.

The equilibrium between the electric forces, the gradient of elastic forces, and the friction forces of chain segments in the gel pore is described by the nonlinear diffusion-type equation.

The "hernias" play the special role in the chain motion. The competition of hernias permits to explain the noticeable orientation of a chain in the field direction, even in rather weak fields. As extensive calculations showed, the deep minimum in "antiresonance" can be obtained only if one takes into account the elastic forces and the non-local interaction of chain segments with one another in chain hernias.

We believe that this theoretical and calculation results are rather important for the understanding of the motion of long DNA chains in pulsed field electrophoresis. They permitted to draw a more accurate physical picture of the motion of the charged polymer chain in an electric field and revealed the special role the hernias and internal friction play for this motion.

The comparison of calculations with experiment showed that the suggested model permits to explain and to quantitatively interpret such anomalies of the DNA mobility in FIGE as "antiresonance" and band inversion, and also to explain the behavior of instant velocities during the periods of positive and negative pulses of an electrical field.

The extensive numerical calculations proved that the suggested "hydrodynamic" model of the polymer chain motion permits to qualitatively and quantitatively describe experiments with gel-electrophoresis and can be used for a prediction of mobility minimums and areas of band inversion, thus enabling us to develop new methods of separation of long DNA chains.

Acknowledgements

We would like to thank Professors G.L. Kotkin, J.-P. Righetti and I.V. Khmelinski for valuable advice and fruitful discussions and support.

References

[1] Schwartz D.C., Cantor C.R. Separation of yeast chromosome-sized DNAs by pulsed field gradient gel electrophoresis // Cell. 1984. Vol. 37. P. 67-75.

[2] Carle G.F., Frank M., Olson M.V. Electrophoretic separation of large DNA molecules by periodic inversion of the electric field // Science. 1986. Vol. 232. P. 65-68.

[3] Sabanayagam C.R., Holzwarth G. Real-time velocity of DNA bands during field-inversion gel electrophoresis // Electrophoresis. 1996. Vol. 17. P. 1052-1063.

[4] Frumin L.L., Peltek S.E., Zhilberstein G.V. Nonlinear focusing of DNA macromolecules // Phys. Review E. 2001. Vol. 64. P. 021902.1-021902.5.

[5] Viovy J.-L. Electrophoresis of DNA and other polyelecrtolytes: Physical mechanisms // Rev. Mod. Phys. 2000. Vol. 72, N 3. P. 813-872.

[6] Deutch J.M. Theoretical studies of DNA during gel electrophoresis // Science. 1988. Vol. 240. P. 922-924.

[7] Lumpkin O.J., DeJardin J., Zimm B.H. Theory of gel electrophoresis of DNA // Biopolymers. 1985. Vol. 24. P. 1573-1593.

[8] Slater G.W., Noolandi J. New biased-reptation model for charged polymers // Phys. Rev. Lett. 1985. Vol. 55. P. 1579-1582.

[9] Duke T.A.J., Semenov A.N., Viovy J.-L. Mobility of a reptating polymer // Phys. Rev. Lett. 1992. Vol. 69. P. 3260-3263.

[10] Duke T.A.J., Viovy J.-L., Semenov A.N. Electrophoretic mobility of DNA in gels I: New biased reptation theory including fluctuations // Biopolymers. 1994. Vol. 34. P. 239-248.

[11] Duke T.A.J., Viovy J.-L. Simulation of megabase DNA undergoing gel electrophoresis // Phys. Rev. Lett. 1992. Vol. 68. P. 542-545.

[12] Duke T.A.J., Viovy J.-L. Motion of megabase deoxyribonucleic acid during field-inversion gel electrophoresis: Investigation by nonlocal Monte Carlo //J. Chem. Phys. 1992. Vol. 96, N 11. P. 8552-8563.

[13] Barkema G.T., Marko J.F., Widom B. Electrophoresis of charged polymesr: Simulation and scaling in a lattice model of reptation // Phys. Rev. E. 1994. Vol. 49, N 6. P. 5303-5309.

[14] Zimm B.H. Size fluctuations can explain anomalous mobility in field-inversion electrophoresis of DNA // Phys. Rev. Lett. 1988. Vol. 61. P. 2965-2968.

[15] Zimm B.H. "Lakes and Strates" model of field-inversion gel electrophoresis of DNA //J. Chem. Phys. 1991. Vol. 94, N 3. P. 2187-2206.

[16] Lim H.A., Slater G.W., Noolandi J. A model of the DNA transient orientation overshoot during gel electrophoresis // J. Chem. Phys. 1990. Vol. 92, N 1. P. 709-721.

[17] DE Gennes P.G. Reptation of a polymer chain in the presence of fixed obstacles //J. Chem. Phys. 1971. Vol. 55, N 1. P. 572-579.

[18] Doi M., Edwards S.F. The Theory of Polymer Dynamics. Oxford: Oxford Univ. Press, 1986.

[19] Grosberg A.Yu., Khokhloy A.R. Statistical Physics of Macromolecules. M.: Nauka, 1989. 370 p. (English translation: 1994, Amer. Institute of Phys.).

[20] Landau L.D., Lifshitz E.M. Statistical Physics. Pt 1. M.: Nauka, 1995. (English translation of a previous Edition:Pergamon, Oxford, 1980) Pt 1. Pergamon, Oxford, 1980.

[21] Riohtmyer R.D., Morton K.W. Difference Methods for Initial-Value Problems, Interscience. 2nd edition. N.Y., 1967.

Received for publication october 25, 2004