Physica A 194 (1993) NorthHolland
163172
Reptation
models for electrophoresis
A. Kooiman
and J.M.J.
InstituutLorentz,
Nieuwsteeg
van Leeuwen
18, 2311 SB Leiden,
The Netherlands
Recently constructed models for reptation such as the RubinsteinDuke model are discussed in relation to electrophoresis. A full solution for the probability distribution for the steady state is given for the case of periodic boundary conditions. The model is modified by allowing only for single and double occupancy of the cells and for this case the drift velocity (diffusion coefficient) is calculated in lowest order of the driving electric field. The influence of the boundary condition is discussed.
1. Introduction In electrophoresis experiments of, e.g., DNA molecules one uses the fact that the mobility of these long polymers is sensitively dependent on their length. The polymers are entangled in a gel (agarose) and move through the gel by reptation. This is a highly complex motion, introduced by de Gennes [l], which is characterised by the fact that the surrounding acts as a fixed confining tube and that the polymer chain first has to store an amount of polymer length before this diffuses to a subsequent position along the tube. The typical experimental situation is that of polymers with a large number of segments N which are driven through the gel by a weak electric field E. So it is important to know how the drift velocity depends on N and E. Particularly interesting and experimentally relevant is the socalled scaling limit in which N+ cfi and E + 0 such that E*N stays finite [2]. Rubinstein and Duke have modelled the reptation by discretizing the motion [3,4]. They take the gel pores as the cells of a regular lattice. The polymer is seen as a chain of segments, to be called reptons, of the size of the persistence length, such that the reptons can be considered as independent units of motion. The reptons hop stochastically from cell to cell, one at every time step. The connectivity of the chain requires that two adjacent reptons are either in the same or in neighboring cells. Initially the string occupies a random sequence of cells (a channel). The internal reptons diffuse in this channel without changing it; in the version of Duke the reptons at the end of the chain may retreat to the inward cell, thereby shortening the channel, or enter a new cell and make the channel longer, all as far as the connectivity of the chain permits. 03784371/93/$06.00
0
1993  Elsevier
Science
Publishers
B.V. All rights
reserved
164
A. Kooiman,
J.M.J.
van Leeuwen
I Reptation models for electrophoresis
We will discuss the RubinsteinDuke model by focussing on the motion of the internal reptons. This is achieved by imposing periodic boundary conditions, turning the endpoints into internal reptons. Thus a substantial simplification of the master equation for the motion is obtained. In fact it allows for a complete solution of the steady state probabilities from which the drift velocity and the diffusion constant can be calculated. The model, as originally proposed, shows an unphysical feature by allowing an unlimited number of reptons to occupy the same cell. The remedy is to limit this number by further restricting the motion. Unfortunately the complete solution for the steady state cannot be found anymore, but an expansion of the drift velocity in powers of the driving electric field E is still possible. A major point of concern is of course the influence of the endpoints. In contrast to the usual thermodynamic limit, boundary conditions are important in the limit of Nco. We will give evidence that in lowest order in the field E the endpoint effects are small, probably negligible in the limit N+ ~0. They are, however, important in the scaling limit.
2. The model In fig. 1 we have drawn a chain on a d = 2 lattice. The direction of the field is along a cell diagonal, such that of two adjacent cells one has a lower potential than the other. The only dynamical entries are the coordinates of the reptons along the field direction. Thus the chain is mapped on the sequence of ldimensional cell coordinates x1, . . . , xN of the cells in which the N reptons are located. The connectivity requires that (xi  x, + , ) is either a, 0 or a (a is the lattice distance). So it is convenient to take as specification of the chain the
E
1‘
Fig. 1. Typical configuration of the chain of reptons. given by y = (0, 1, 1, 1, 1, 0, 1, 1, 1.0, 1).
The internal
coordinates
of the state
are
A. Kooiman,
J. M. .I. van Leeuwen
I Reptation models for electrophoresis
165
center of mass coordinate
(2.1) and the differences YiCxi+l
 x,)/a
)
(2.2)
which are restricted to the values  1, 0, 1. Translational invariance in space permits us to discuss the motion of the center of mass by Fourier decomposition, such that we can concentrate for the spatial homogeneous case on the motion of the states y = y,, . . . , y,_, . Transitions which are allowed by the connectivity have a nonzero transition rate
where w is a hopping rate and B a bias B = exp(aEq/2k,T)
= exp(e/2),
(2.4)
with T the temperature, q the charge and k, the Boltzmann constant. Only one repton hops at a time and W+ applies when the repton moves in the direction of the field and W_ when it moves against the field. The probability P( y, t) of state y develops in time according to the master equation
r)
WY,
=
dt
c [WYIY’)
P(Y’Y
ys
The steady state distribution
F [WYIY')
JYY')

9
WY'IY)
P(y) therefore
WY'IY)
P(Y)1
P(Y?
[)I
’
(2.5)
obeys
=0 *
(2.6)
Once P(y) is known the drift velocity u can be calculated as (2.7)
where u(y) is the velocity in the state y computed U(Y) =
cY' W+(Y’lY)

w(Y’lY)l
7
as
(2.8)
166
A. Kooiman,
J. M.J.
van Leeuwen
I Reptation models for electrophoresis
being the difference of the possible up and down moves. The computation (2.8) of u(y) is simple, the determination of P(y) from (2.6) is the hard part of the problem and the summation (2.7) is still a substantial task.
3. Periodic boundary conditions In order to introduce the periodic boundary conditions we add one extra variable yN = 1, 0, 1 to the chain, specifying the position xN+, of repton N + 1 with respect to N. The motion of repton N + 1 is simultaneous and in the same direction as repton 1. One can see the chain as repeating itself with N + 1 the image of 1. x,,,+,  x, = a C, y, gives the difference in position of the images along the field direction. Thus all yi become equivalent and, more important, the whole channel which the chain traces out becomes an invariant. Going along the chain from y1 to y,, every yi # 0 means a connection between two adjacent cells. The values yi = 0 represent the stored length: a single yi = 0 means a doubly occupied cell (by repton i and i + l), a pair y, = yj + 1 = 0 means a triply occupied cell (by i, i + 1 and i + 2), etc. The stored length (the zeroes in y) diffuses through the channel without altering it. So the collection of yi # 0 does not change during the motion, only their indices change. Therefore one can take a different representation of the state y by considering the collection s, , . . , sL with si = ?l for the set of nonzero y’s. S, positions cell 2 with respect to 1 and so on while sL tells the position of the periodic image of cell 1 with respect to cell L. The zeroes among the y’s are represented by an occupation n( = 1, 2, . . , which counts the number of extra reptons (extrons) in cell 1. The connectivity requires always at least one repton in a cell. The motion is a redistribution of the N  L extrons. Consequently the master equation can be reformulated in terms of the probability P(n) for the occupation n,, . . . , nl_ given a channel sir . . . , s,. The steady state equation reads
~~
[e(ni) Bsi + 8(n;+l) B“‘]‘(n,,
..
1
nL>
= g, [@z,) BmSfP(n,, . . , ni  1, IE;+, + 1, + O(n,+,)
BS’P(n,,
. , ni +
. . . , nt_)
1, Al;+,  1, . . . , n,)] .
(3.1)
The lefthand side gives the gain terms when extrons move between cell i and i + 1, whereas the righthand side gives the corresponding loss terms for the probability. The step function 0(n),
A. Kooiman,
O(n)=
J. M. J. van Leeuwen
0,
n=O,
1
na1,
1?
I Repration
models for electrophoresis
167
(3.2)
has to be inserted because the processes can only take place when n > 0. It is instructive to inspect (3.1) for the case that only one extron is available (L = N  1). Let pi be the probability that the extron is in cell i. The master equation reduces in this case to Pig
s~l
_Pi_,B”zl
=Pi+,BSc
piB”t
,
(3.3)
which is a simple recursion relation with periodic boundary condition pL+ 1 = p,. It can be solved by a technique due to Derrida [5], which gives the pi as explicit functions of the channel parameters s, , . . . , sL. Now the remarkable property of (3.1) is that the solution of P(n) can be expressed in terms of the pi as [6] (3.4) Viewing p, as a Boltzmann factor pt = exp( PE,), one may interpret (3.4) as the distribution of independent bosons occupying the levels E,. The factor +(s) gives the amplitude of the channel s. For our periodic boundary conditions the channels s are invariant and $(s) is a free quantity. Assuming equal probabilities for all channels, the drift velocity has been calculated and a scaling form results of the type (3.5) where z = NE’. The scaling function Q(z) is plotted in fig. 2. The values G(O) = 1 and Q’(O) =  & have been determined analytically [6]. By finite size scaling G(z) has been determined accurately up to arguments z = 40, using chains of length up to N = 30 [7]. The full curve is an approximation, discussed in ref. [6].
4. Fermionic extrons As we pointed out, the extrons behave as independent bosons, having a tendency to condense in the lowest level, i.e. in the cell with the highest single extron probability pt. Such crowding of the cells, which becomes stronger with increasing fields, is clearly unphysical, given the fact that the gel pores are of
168
A. Kooiman,
.I. M. J. van Leeuwen
/ Reptation models for electrophoresis
I
x
z
.2 
0 0
I
I
I
I
20
40
60
60
100
*
Fig. 2. The scaling function Q(z) for the RubinsteinDuke condition and equal probabilities for the channels. The points scaling and the full curve is an approximation.
model with periodic boundary have been obtained by finite size
similar size as the reptons. Thus we propose to further limit the motion by restricting the occupancy of the cell. Most natural and simple is to allow the cells to be occupied by zero or one extron (one or two reptons), which implies that the extrons behave as fermions. The equivalent of (3.1) for the fermions can be written as
,$,[n,(l  n;+l)BsI + n;+,(l=i
[n,(l
n;P"'lJYn,,~
ni+,)B“’ + n,+i(l
.
2
nt*1
n;)W] (4.1)
The factors n,(l  IZ,+~) and (1  lzi)ni+, enforce the occupation to obey the exclusion principle. The 1extron problem is of course the same as (3.3) but the product property does not apply to the many extron probability. In spite of the fact that an explicit expression for P(n) is lacking, we can make an expansion for the drift velocity in powers of E, by taking moments of the master equation (4.1). Summing (4.1) over all n leads to the identity, expressing the conservation of probability. Multiplying (4.1) with factors nk leads to moment equations after summation over n. The simplest result follows from a single factor nk, leading to the relation [7]
A. Kooiman,
J. M.J.
van Leeuwen
(JkI(%4) = (Jkh 4)
I Reptation models for electrophoresis
169
(4.2)
7
where Jk(n, s) is given by J,(n, s) = w[n,(l
nk+,)BS*  (1
n,)n,+,BP]
.
(4.3)
( ) stands for an average with P(n) as weight, Jk is the net flow from cell k to k + 1 and thus (4.2) expresses the fact that this flow
(Jkh s)) = J(s)
(4.4)
is constant along the channel. With this interpretation average total velocity in state (n, s) as
one can express the
where S = C, sk is the endtoend distance of the polymer along the direction of the electric field. Now (4.4) is powerful enough to determine the lowest order of J(s) in the field. Expanding B = 1 + e/2 in (4.3) yields for (4.4) (4.6) Summing (4.6) over k, the first term drops out and in the second term one may compute the average to lowest order in the field (B = l),
bk +
nk+l
consequently

2nknk+lh
=
2(N  L)(2L  N) . ’ L(L  1)
(4.7)
(u(n, s)) can be found to lowest order in E as [7]
(4% 4) = w
w

J5)W

w
s2E
L2(L1)
The sum over channels s can be computed, the channels, with the relation T s2 = L2L.
(4.8)
.
assuming equal probabilities
for all
(4.9)
Then, by doing the sum over L by steepest descent, the lowest order for the drift velocity becomes u = (awlN)(3fl
5)e
(4.10)
170
A. Kooirnan, J.M.J.
van Leeuwen
1 Reptation models for electrophoresis
for large N. Comparing this to (3.5) we see that the coefficient of the linear term in E, which is f for the bosons, is replaced by 3fl5 = 0.196 for the fermion case. Not surprisingly the drift is reduced by the restriction on the motion by the fermion exclusion principle. Higher orders in E may be obtained by taking higher moments of the master equation (4.1) but the procedure becomes quite involved as the next order is already E’. The moment method works also for bosons, but it is essential that one can take moments at fixed channels s. The moment method applied to the general equation (2.6) yields a number of interesting relations, but not an explicit expression for the lowest order drift velocity, which is by Einstein’s relation equivalent to the diffusion coefficient in zero field.
5. Effects of boundary
conditions
An important question is the influence of the boundary condition on the motion of the polymer chain. Let us compare chains with an equal number of variables y,, . . . , y,, one being open and having N + 1 reptons and the other being periodic with repton N + 1 moving as image of repton 1. It is clear that the latter can only perform a subset of the motions of the open chain. On the basis of the principle that more freedom will allow a faster motion one has
(u/v+ I Jopen s
(UN)periwiic
’
(5.1)
The missing moves of the periodic chain cause the ystate space to break up into a number of subspaces, which do not communicate with each other. The subspaces correspond to the invariant channels. As mentioned before, the relative importance of the channels is free in the periodic chain. The open chain yields the probability for each state y, so one can calculate the probabilities of the channels as well as the probability distribution inside the channel. We have done so for small chains and conclude that the distribution in a channel does agree well with that following from the periodic chain, at least for weak fields. The expression (4.8) shows that to lowest order in the field, only the zero field channel probabilities $(s) enter in the expression for the drift velocity. So if the channels give the proper contribution to the drift, the periodic chain will give for E + 0 the correct expression (3.5) and (4.10). In physical terms it means that the drift is limited by the internal reptons and not speeded up for large N by the larger freedom of the endpoint motions. In fig. 3 we have plotted the scaled diffusion constant for the fermionic reptons as it can be evaluated for short chains (up to N = 12) by direct solution of the master equation. Since the steady state is given by the eigenvector
A. Kooiman,
J.M.J.
van Leeuwen
I Reptation
models for electrophoresis
171
N
Fig. 3. The scaled diffusion coefficient versus the number of reptons for the fermion case with free endpoint motion. The values have been determined by a direct solution of the master equation for chains up to N = 12 and by a Monte Carlo simulation for larger systems. The error bars give the standard deviations and the dotted line is the asymptotic value N*Dla’w = 3fi  5 for the periodic chain.
corresponding to the largest eigenvalue it can be obtained by repeatedly applying the matrix of the master equation to an arbitrary vector. For larger N a few points have been determined by Monte Carlo simulation. One observes that the values approach the asymptotic value 3V?!  5 of the periodic chain. For a definite conclusion one would need more accurate large N values, which are difficult to obtain. In this respect the situation is similar to the bosonic case
PI. For stronger fields, effects of orientation start to play a role. Each time a new cell is probability B si by the first repton and with So an estimate for the channel probabilities
h,(s) =
(;)(B”+B?,
of the polymers by the field will added to the chain it is done with probability B”l by the last repton. is (5.21
172
A. Kooiman,
J. M. J. van Leeuwen I Reptation models for electrophoresis
where S = C, si, showing that the more oriented channels acquire a larger weight. Using this ansatz we can repeat the calculation for the drift velocity, and for the boson chain we obtain (5.3) Note that this has the opposite sign in third order in E in contrast to (3.5). Therefore the drift velocity will level off for longer chains and becomes less sensitively dependent on the length of the polymer, a phenomenon observed in experiments. We find a similar behavior for the fermionic open chain from exact solutions for short chains. Widom et al. predicted that the scaling function will be dominated for large fields by this third order term. For large fields one finds however that the drift velocity vanishes exponentially such that the scaling function, after going through a maximum, decays to zero. A fuller account will be given elsewhere [7].
References [II P.G. de Gennes, Scaling Concepts in Polymer Physics (Cornell Univ. Press, Ithaca, Chem. Phys. 55 (1971) 572. J. Phys. I (Paris) 1 (1991) 1759. 121 B. Widom, J.L. Viovy and A.D. Defontaines, Phys. Rev. Lett. 59 (1987) 1946. 131 M. Rubinstein, 141 T.A.J. Duke, Phys. Rev. Lett. 62 (1989) 2877. J. Stat. Phys. 31 (1983) 433. 151 B. Derrida, Physica A 184 (1992) 79. 161 J.M.J. van Leeuwen and A. Kooiman, J. Chem. Phys., to be published. 171 A. Kooiman and J.M.J. van Leeuwen,
1979); J.