Collapse of a semiflexible polymer in poor solvent
Alberto Montesi,1,2 Matteo Pasquali,1,2,* and F. C. MacKintosh2,3
1Department of Chemical Engineering, Rice University, 6100 Main Street, Houston, Texas 77005, USA2The Kavli Institute for Theoretical Physics, University of California, Santa Barbara, California 93106, USA
3Division of Physics and Astronomy, Vrije Universiteit, 1081 HV Amsterdam, The Netherlands~Received 29 August 2003; revised 3 December 2003; published 27 February 2004!
We investigate the dynamics and pathways of the collapse of a single, semiflexible polymer in a poor solventvia three-dimensional Brownian Dynamics simulations. An example of this phenomenon is DNA condensationinduced by multivalent cations. Earlier work indicates that the collapse of semiflexible polymers genericallyproceeds via a cascade through metastable racquet-shaped, long-lived intermediates towards the stable torusstate. We investigate the rate of decay of uncollapsed states, analyze the preferential pathways of condensation,and describe the likelihood and lifespan of the different metastable states. The simulations are performed witha bead-stiff spring model with excluded volume interaction, bending stiffness, and exponentially decayingattractive potential. The semiflexible chain collapse is studied as a function of the three relevant length scalesof the phenomenon, i.e., the total chain lengthL, the persistence lengthLp , and the condensation lengthL0
5AkBTLp /u0, whereu0 is a measure of the attractive potential per unit length. Two dimensionless ratios,L/Lp andL0 /Lp , suffice to describe the dimensionless decay rate of uncollapsed states, which appears to scaleas (L/Lp)1/3(L0 /Lp). The condensation sequence is described in terms of the time series of the well separatedenergy levels associated with each metastable collapsed state. The collapsed states are described quantitativelythrough the spatial correlation of tangent vectors along the chain. We also compare the results obtained with alocally inextensible bead-rod chain and with a phantom bead-spring chain. Finally, we show preliminary resultson how steady shear flow influences the kinetics of collapse.
The conformation of individual polymer chains depenon the solvent properties@1–3#. In good solvent, the monomers effectively repel each other, preferring to be surrounby solvent. This effect leads to a swollen coil conformatifor flexible polymers in good solvent. In poor solvent, coversely, the monomers try to exclude the solvent and eftively attract one another, and a flexible chain forms a copact globule of roughly spherical shape to minimize tcontacts between monomers and solvent. The dynamicthe coil-globule transition in flexible chains are well know@4–7# and involve the formation of a pearl necklace andgradual diffusion of large pearls from the chain ends@8,9#.
An early theoretical picture of the phenomenon predicthe initial formation of a dense, uniform, sausage-like shaggregate, which in a second stage forms a globule, drby surface area minimization@4#. The first experimental evi-dence of a two-stage kinetic was later obtained using acapillary tube cell for dynamic light scattering measureme@8#. More recently, simulations showed that capillary insbility prevents the formation of such sausage-like shapesthat the condensation pathway actually involves an inifast crumpling of the unknotted polymer chain, with the fomation of pearls of collapsed, dense phase connected bycollapsed bridges, and a subsequent slow rearrangemethese pearls to form a compact globule@6,9#. Each stage ischaracterized by a length scale, and by a time scale incring with the relevant length scale. The pearling stage
volves only local rearrangement of the chain configuratiowhile the dimensions and configurations of the whole chare weakly modified. The pearls then grow attracting momers from the bridges, which become increasingly stretcand then shrink as the collapse proceeds. Eventually,pearls come into contact, coalescing into a single globule
Many polymers, however, exhibit substantial bendistiffness, i.e., they are semiflexible, and they collapse in psolvent towards very different equilibrium states, and folloing different kinetics. Semiflexible polymers can be dscribed by the wormlike chain model; examples inclubiopolymers~such as F-actin, DNA, and xanthan! as well assynthetic polymers~such as PPTA and PBO! used for pro-duction of high-performance fibers~e.g., KevlarTM andZylonTM). Such polymers form open, extended structuresgood solvents. In poor solvent, the effective self-affinitythe chain leads to chain collapse. The collapsed statefigurations and the pathways to their formation are the reof the interplay between two opposing forces: the bendforce related to the chain stiffness and the attractive fodue to the poor solvency of the environment.
A compact globule is energetically unfavorable for semflexible polymers because it involves bending over shlength scales. The collapsed ground state is instead a towhich reduces the monomer-solvent contacts without caing excessive bending penalty, i.e., ensuring that the chastill roughtly straight on short length scales. Theoreticanalysis@10,11#, Monte Carlo simulations@12#, and experi-mental evidence@13,14# have detailed the stability, the features, and the packing of the collapsed tori. The dynamicthe collapse of wormlike chains in poor solvent have beinvestigated more recently@15–17#, strongly suggesting a
MONTESI, PASQUALI, AND MACKINTOSH PHYSICAL REVIEW E69, 021916 ~2004!
possible generic pathway for the collapse of semiflexipolymers, featuring a series of long-lived, partially colapsed, racquet-shaped intermediate states, before thetual collapse to a torus. Such intermediate states formenergetically driven cascade of increasingly compact conmations with sharp transitions between their energy leve
The dynamics of collapse of semiflexible moleculesrelevant for understanding DNA condensation. DNA frquently forms condensed structuresin vivo, e.g., in DNAreplication, viral transfection, and compaction within speheads and nucleosomes@14#. Controllablein vitro DNA con-densation is involved in nonviral gene therapy—whichused in the treatment of various kinds of disease, includcancer—for packaging in, and transfection from gene deery vehicles. Multivalent cations, in particular polyaminehave been used to induce condensationin vitro, showing thatcondensed DNA exists in toroidal and rodlike structurPolyamines present ubiquitously in living cells in millimolaconcentrations~e.g., putrescine and spermidine! are also be-lieved to induce condensationin vivo @18#. A better under-standing of the dynamics of condensation would lead toprovements in gene delivery technology, as it may hdefining the best conditions and protocols to obtain DNcondensation.
In the present work, we verify previous hypotheses@15#of the suggested pathway for the collapse of semiflexpolymers and further explore the effect of the relevant lenscales on the phenomenon. The natural length scale fopolymer stiffness is its intrinsic persistence lengthLp , de-fined as the ratio of the bending stiffnessk to the thermalenergykBT, which represents the length along the chainwhich tangent vectors remain correlated in a theta solvi.e., in absence of energetic~attractive or repulsive! interac-tions between chain and solvent. In a poor solvent, theance between intrachain attractive forces and bending focan be considered in terms of the so-called ‘‘condensalength’’ L05AkTLp /u0, whereu0 is the value of the attractive potential per unit length. We choose as time scale forcollapse event the rotational diffusion time of a perfecrigid rod with L5Lp , t rot,Lp
5zLp3/(72kBT), wherez is the
transverse friction coefficient@19#. Therefore we uset rot,Lp
as unit of time andkBT as unit of energy in presenting all thresults of the present work. We investigate the sequencesthe likelihood of different collapsed structures formedvarious solvent quality~described byu0) and chain stiffness
II. SIMULATION DETAILS
The dynamics of the continuous wormlike chain are dscribed by an inertialess Langevin equation:
z•H ]r
]t2g•r J 52
dU
dr1h. ~1!
Equation~1! balances the hydrodynamic force exerted onchain by the surrounding fluid with the force due to tglobal potential acting on the chainU5Ubend1Uconn1Usolv , and the Brownian force per unit lengthh, withcorrelations ^h(s,t)h(s8,t8)&52kBTzd(t2t8)d(s2s8).
02191
e
en-n
r-.
-
g-,
.
-p
ehhe
rt,
l-esn
e
ndt
-
e
Ubend5kTLp*0Lds]u(s)/]s is the bending potential relate
to the chain stiffness,Uconn is the connector potential,Usolvdescribes the monomer-monomer interaction, andBrownian term accounts for the rate of fluctuating mometum transfer due to the random collisions of small solveparticles with the chain.
The potentialUsolv chosen to account for the effect opoor solvent is an exponential attractive potential withexcluded volume repulsive term:
Usolv52E dsE ds8u0expS 2ur ~s!2r ~s8!u
RattrD
1E dsE ds8u0S s
ur ~s!2r ~s8!uD 12
, ~2!
whereu0 defines the depth of the potential well,Rattr is therange of the attractive forces, ands is the radius of theexcluded volume region. With this potential,u05Rattru0.Different choices for the shape of the potential shouldchange the stable and metastable configurations forwormlike chain in poor solvent, whereas they modify tactual value of the condensation lengthL0 and they may alterthe speed and the dynamics of the collapse. Our choice oexponential potential is related to the physics of the condsation phenomenon. DNA and polyelectrolytes collapsecause of charge inversion and counterion induced attrac@20#. The attractive forces are therefore due to highscreened electrostatic charges, which typically show anponential decay, mimicked by the potential selected in twork.
The Brownian Dynamics code used in this work solvesdimensionless, discretized form the Langevin equationmotion of a linear wormlike chain ofN beads with positionsR1 , . . . ,RN , connected byN21 connectors of equilibriumlength a[L/(N21) with unit tangent vectorsui[(Ri 112Ri)/a. For most of the simulations reported in this worwe use quadratic springs as connectors, with a force consequal to 100kT/a2, which ensures a small variation of thconnectors’ length. This choice allows for longer time stewithout interfering with the spring relaxation time. We usemidstep algorithm@21# to compute the bead trajectories geerated by the equation of motion
zbH dRi
dt2g•Ri J 5Fi5Fi
bend1Fielast1Fi
solv1Firand . ~3!
Here,zb5za is a bead friction coefficient, the bending forcis
Fibend5k/a (
j 52
N21
]~uj•uj 21!/]Ri , ~4!
the elastic force due to the spring is
Fielast5H@~ uRi2Ri 21u2a!ui 21
2~ uRi 112Ri u2a!ui #, ~5!
6-2
beor
aFnaths
adfo
d
haf
deee
-phethiseto
senethayt thu
rtsheicrinceatnbt
froednis
esrent
heII E.
mi-ol-ientu-
ni-eti-ith
thepseob-o
ectheare
iony totu-an-no-ly,ateichn,ing.o-
ng
COLLAPSE OF A SEMIFLEXIBLE POLYMER IN POOR . . . PHYSICAL REVIEW E69, 021916 ~2004!
the potential force related to the effective interactionstween the beads of the chain, repulsive within the hard cattractive at short distance is
Fisolv5u0(
j Þ iF expS 2
uDRi j uRattr
DRattr
212s12
uDRi j u13G DRi j
uDRi j u, ~6!
whereDRi j 5Ri2Rj , and the random force isFirand .
Another relevant issue in the simulations is the accurdescription of the high curvature of the collapsed states.a fixed level of arclength discretization, the computatioerror in the discretized wormlike chain model grows aschain curvature increases. To achieve an accurate repretation of highly curved states characteristic of raquet heand tori ~see Sec. III below! a sufficiently large number oconnectors within the characteristic length scale of the clapsed structure is needed. The typical size of the tori anthe racquet heads are of the same order of magnitude ofLp .In the simulation conditions examined here, we found tusing 20 connectors perLp yields an accurate description ochain curvature.
We also performed simulations with a bead-rod mo@22,23#, which ensures local inextensibility of the wormlikchain. The computational cost for this model is, howevsignificantly higher (;10–50 times! than that for a beadspring chain model, mainly because a smaller time sterequired for convergence. While performing most of tsimulations using the bead-spring model, we comparedresults for few conditions with the bead-rod model, as dcussed in Sec. III D. In that section we also compare somthe results with that obtained using a bead-spring phanchain, i.e., without excluded volume interactions.
For simplicity, we use isotropic drag in the simulationand we neglect hydrodynamic interactions between segmof the chain. Including such interactions would likely givminor corrections to the evaluation of the time scales ofphenomena, but it would not change the kinetic pathwand the lifetimes of the intermediate states, because abeginning of the simulations the chains are locally straigand after the collapse the hydrodynamic interactions wobe dominated by the attractive forces.
III. RESULTS
This section is organized as follows: Sec. III A repoqualitatively the results of our simulations, describing tvarious possible intermediate configurations and the typpathways. Section III B presents a more quantitative desction of the collapsed states in terms of the correlation fution between tangent vectors along the chain and of theergy levels associated with each different collapsed stwhich confirm that the toroid is the stable configuration, ashow how the collapse proceeds through a cascade of suquently more energetically favorable metastable states,various multiple-heads racquets shapes, as predictedtheory @24#. The decay rate of the number of uncollapsstates is investigated as a function of a proper combinatioL, L0, andLp in Sec. III C, together with a statistical analys
02191
-e,
teorleen-s
l-of
t
l
r,
is
e-ofm
,ts
eshet,ld
alp--n-e,dse-hem
of
of the preferred collapsed configurations. Sec. III D provida comparison between the results obtained with the diffecomputational models~phantom chain and bead-rod chain!.Some preliminary results of the effect of shear flow on tkinetics and pathways of collapse are presented in Sec. I
A. Molecule configurations and collapse pathways
Figure 1 shows a typical snapshot of an ensemble of seflexible chains during collapse. The 3D structures of the mecules are shown here as 2D projections on a convenplane. At a given time molecules coexist in different configrations, from completely uncollapsed@molecule~a!#, to par-tially collapsed@molecules~b!, ~f!, and~g!#, to higher orderracquet@molecules~c! and ~e!#, to torus@molecule~d!#.
The fraction of collapsed molecules increases monotocally with time, confirming that collapsed states are energcally favorable in a bad solvent. Many other ensembles wdifferent values ofL/Lp andL0 /Lp show a similar behavior.However, if the attractive forces are not strong enough orbending stiffness is too high there is no evidence of collaand the open conformations persist in time. As alreadyserved in@15#, a collapse sequence can evolve following twdifferent generic patterns. The first pattern involves the dirformation of a torus from the open conformation, when ttwo ends of the chain meet while the tangent vectorslocally parallel and pointing in the same direction (u1•uN'1). In this instance, the chain collapses without formatof any intermediate, the end to end distance drops quickla value of the order of the torus dimension, and then flucates as the chain keeps folding into the final structure. Mewhile the total energy of the chain keeps decreasing motonically until it reaches the equilibrium value. Conversethe second pattern involves the formation of intermedimetastable states: initially a single-headed racquet whfolds quickly into a short-lived, 3-head racquet configuratioand then rearranges as a toroid via a subsequent foldFigure 2 shows the two different pathways. We n
FIG. 1. Typical configurations of an ensemble of collapsichains at an intermediate time (L/Lp58, u051.25kT/a).
6-3
uncetaw
ode
o
nin
blre
ther-
oseles
en-ev-wellhest
ula-h.s
rmalcaseticeen
thetive
in isgu-odrac-
lev-pan-ay
u-n-ntforbili-u-
t:ionbe-the
tod
MONTESI, PASQUALI, AND MACKINTOSH PHYSICAL REVIEW E69, 021916 ~2004!
tice here that the direct formation of a torus is a morelikely pattern in 3D—and in fact observed with scarfrequency—than in 2D. Whereas in 2D the relative oriention of the two ends is described by a single angle, tangles are required in 3D and both need to be close to 2p forthe direct formation of a torus to take place; therefore, mof the chains initially form a single headed racquet, asscribed below@25#.
The time evolution of the total energy of an ensemblecollapsing chains is shown in Fig. 3. The total energyU ofthe chain during collapse is defined as the sum of the being energyUbend and the interaction energy along the chaUsolv . In the discretized form,
U5kTLp (i 51
N22
cos2u i1(j 51
N
(k5 j 11
N
usolv, jk , ~7!
where u i is the bond angle andusolv, jk is the interactionenergy evaluated between the beadsj andk along the chain.
The time series in Fig. 3 reports most of the possievolutions observed in our simulations. Few molecules
FIG. 2. ~Color online! Comparison between time series of endend distanceR and total energy for direct formation of toroid anformation through intermediate states.
02191
-
-o
st-
f
d-
e-
main in the uncollapsed states for the whole duration ofrun, equivalent to 250t rot,Lp
. For these molecules, the the
mal fluctuations never lead two branches of the chain clenough together to initiate the collapse. The other molecutypically reach lower energy states by collapsing in sequtially higher order metastable conformations. The energy lels for each racquet-head shape of increasing order areseparated, and there is a rather large gap from the higorder racquet observed~sixth order! and the stable tori, asexpected from theory@24#. It is interesting to compare theactual values of the energy levels obtained in these simtions with those calculated using the theoretical approac
Schnurret al. @24# obtained the conformational energieof the intermediate racquet states in the absence of thefluctuations, relating the annealed shapes to a particularof Euler’s elastica. In their calculations, the relative energeadvantage for segments of the chain to bundle has bevaluated expressing the energy as a surface tensionG ofsolvent-exposed sites. For the simulation condition ofdata reported in Fig. 3 we get the two sequences of relaenergies~theoretical and from simulations! reported in TableI. In both cases, the energy of a straight, uncollapsed chaset to 1, and we report the energy of the collapsed confirations relative to this one. We note that there is a very goagreement for open chain, 1-head, 2-heads, and 3-headsquets, whereas the simulation results show lower energyels for higher order racquets. The reason for these discrecies, substantial for the more compact conformations, mlie in the finite range of the attractive potential in the simlations. While the theoretical calculation using surface tesion accounts only for the local coverage of the filamecross section, the finite range in the simulations allowsnext-nearest neighbor interactions, leading to higher stazation for tighter structures~where each bead has more nmerous next-nearest neighbors!.
Most of the chains initially form a single head racqueafter the first monomer-monomer contact and the formatof the racquet head, the two sides of the chain tend tocome parallel and form a neck region, which ensures
ns and
FIG. 3. ~Color online! Time series of total energy for an ensemble of molecules, with labels of the corresponding configuratioactual shape of~a! 1 head,~b! 2 heads,~c! 5 heads racquets, and~d! torus.
6-4
ouoennF
diab
asdmtho
epolia
rtaneor
tioth
u
olizleore
ts
ais
ith
d bybe-gleckost
a
toe
as aunc-
sewe
ce
the
ofs isthe
ndre-le
amval-
iss of-tiveith
ity
ng
fivelyThedingtelyrealy-
ofbedge
tas
COLLAPSE OF A SEMIFLEXIBLE POLYMER IN POOR . . . PHYSICAL REVIEW E69, 021916 ~2004!
maximum number of monomer-monomer contacts withincreasing the bending penalty. The molecules then fagain on themselves, forming higher order shapes. Deping on how the molecule folds back on itself, the next cofiguration can be a 2-heads or a 3-heads racquet shape.ure 3 shows examples of both pathways. Successive folleads to higher order racquets. The transition between esingle metastable conformation is very rapid and driventhe deterministic attractive and bending forces. In contrthe initiating event for each successive collapse is relatethe Brownian motion of the molecule, and therefore the tiinterval between folding events is completely random forsimulations that we have performed. In the simulationslonger chains (L510Lp) we observe combinations of thconformational elements described above, e.g., chainstially folded in a racquet shape at one end while still unclapsed at the other. Those conformations have intermedenergies between the well-defined levels correspondingracquets and tori and quickly disappear in favor of mocompact and ordered structures. In contrast, the metashigh order racquets are extremely long-lived, and the egetic barrier for reaching the equilibrium shape is of theder of severalkBT. While the partial unfolding of a torus ina racquet is never observed, simulations show the formaof a torus via metastable configurations, and this, togewith the lower total energy shown by the torus~see Fig. 3!,confirms once more that the latter is the equilibrium configration.
B. Quantitative description of collapsed states
The configurations and the shapes of semiflexible mecules in bad solvent can be studied through direct visuation of the collapsing chain. While this method is a valuabway to gather information on the kinetics and pathwayscollapse, it is time-consuming and not quantitative. We thefore define a spatial correlation matrixM to analyze in amore systematic way the collapsed shapes. The elemenM are defined asMi j 5ui•uj . With this definition,M is sym-metric; a perfectly straight rod hasMi j [1 ; i , j , as all theconnecting normalized vectors have the same directionorientation. An uncollapsed semiflexible chain with perstence lengthLp will have Mii [1, while the off-diagonaltermsMi j would be smaller and decaying exponentially wthe distance along the chain between the connectorsi and j.
TABLE I. Comparison between the energy level of the mestable condensed states from theory@24# and present simulation(L/Lp58.0, L/L0565.3).
When the chain collapses, the actual shape assumethe molecule can be inferred from the spatial correlationtween connectors. The spatial correlation matrix for a sinhead racquet, for example, will be a block matrix: one blowill show values close to 1, corresponding to the almstraight filament in the neck of the racquet, followed byregion with rapidly decreasing values from 1 to21—theracquet head—and by another block with values close21, indicating the other filament of the neck, with the samdirection but opposite orientation.
The spatial correlation matrix can be simplified intospatial covariogram for the molecule. The covariogram istandard statistical measure of spatial covariance as a ftion of distance@26#, and is defined as
C~h!51/N~h!(j 51
N
(i 51
N
$Zi2m%$Zj2m%, ~8!
where h is the distance~in time or in space! between theobservationsZi andZj , N(h) is the number of observationat a distanceh, andm is the expected average value of thobservation. In the present analysis, for each moleculeconstruct the spatial covariogramCsp assumingZi5ui andm50W . We therefore obtain a scalar function of the distanalong the chains, bounded between11 and 21, whichcontains all the relevant information about the shape ofmolecule.
Figure 4 shows three relevant examples of this typeanalysis: the actual shape of the semiflexible moleculeshown together with the spatial correlation matrix andspatial covariogram. For a uncollapsed chain@Fig. 4~a!#, thematrix shows randomly alternating regions of positive anegative correlation along the length of the chain, corsponding to the typical open configuration of a semiflexibchain in good solvent. The corresponding covariogrshows how the connectors correlation decays from 1 toues around 0 along the length of the chain; the decayexponential at short distances and related to the stiffnesthe chain, i.e.,Csp5exp(2s/Lp). At longer distances, the correlation randomly fluctuates between negative and posivalues. As shown in Fig. 4, a single exponential curve wthe theoretical coefficient21/Lp fits very well the spatialcovariogram of the open chain, thus confirming the validof this quantitative analysis.
For a racquet shape with multiple heads,M has a struc-ture of alternating blocks with strong positive and stronegative correlation, as shown in Fig. 4~b!. The five straightfilaments in the neck of the racquet correspond to theseblocks in the matrix; the size of the blocks is approximatethe same and indicates the length of the racquet neck.covariogram shows four zero-crossings, each corresponto a head of the racquet, and can be fitted quite accurawith a fourth order harmonic; the parameters of this fit arelated to the curvature of the racquet heads. From the ansis of the covariogram only, we can extract:~1! the numberof heads in the racquet, which is equal to the numberzero-crossing of the covariogram. This relationship caneasily understood recalling that the tangent vector at the e
-
6-5
MONTESI, PASQUALI, AND MACKINTOSH PHYSICAL REVIEW E69, 021916 ~2004!
FIG. 4. ~Color online! ~a! Actual shape,~b! correlation matrix, and~c! spatial covariogram~with fitting curve! for ~1! an uncollapsedmolecule,~2! a metastable racquet-head shape, and~3! a stable torus.
eccco
a
m-eththn
ai
n
lepri
ano
ate
the.
-e of
bleagech
si-
nd a
.,
tiveluesa-
e
for
of the head is perpendicular to the tangent along the nand~2! the size of the collapsed molecule, from the distanbetween two subsequent zero-crossings. In the collapsefiguration in Fig. 4~a.2!, for example, the length of the4-heads racquet can be correctly estimated from the covogram to be about 2Lp .
For a torus, the spatial correlation matrix appears copletely different@Fig. 4~c!#, showing a structure with diagonal bands. In fact, the correlation between pairs of conntors at the same distance is the same, regardless ofposition along the chain; this translates into bands alongdiagonals of the matrix. The distance between a diagowith values close to 1 and one with values close to21corresponds to the diameter of the collapsed torus. The sinformation can be extracted from the covariogram, whicha perfect sinusoid, instead of a higher order harmonic fution: the period of the sinusoidal function corresponds toptimes the diameter of the torus, and the number of compperiods in along the chain indicates the number of looformed by the collapsed molecule. Again, from the covaogram alone we can derive~1! the number of loops formedby the torus, equal to half of the number of zero-crossing~2! the size of the collapsed chain, which is in the caseFig. 4~a.3! equal to roughly 1.5Lp .
C. Kinetics of collapse
The fraction of collapsed semiflexible molecules in a bsolvent increases monotonically with time, as already sta
02191
k;en-
ri-
-
c-eire
al
mesc-
tes-
df
dd.
We have studied the effect of the relevant length scales ofphenomenon,Lp , L0, andL, on the kinetics of this collapseBy systematically varying the total length of the chain (L53Lp , 5Lp , 8Lp , and 10Lp), as well as the strength (u050.25kT/a2, 0.5kT/a2, and 1kT/a2) and the range (Rattr55/40Lp and 1/6Lp) of the attractive forces, we have explored the effect of these parameters on the decay ratuncollapsed molecules. We define the time to collapsetcoll asthe time required by the molecule to form the first metastacollapsed structure. We obtain from simulations the avervalue of tcoll in an ensemble of collapsing chain under eacondition, normalized witht rot,Lp
. The inverse oftcoll is
taken as a measure of the decay rateDcoll of uncollapsedmolecules. The distribution of the times to collapse followthe expected lognormal distribution, typical of such dynamcal processes, with a peak close to the average value along tail, as shown in Fig. 5.
At a given value ofL/Lp , the decay rate increases, i.ethe molecules collapse faster, for decreasing values ofL0,i.e., for increasing values ofu0. Similarly, the decay rateincreases for longer chains at a given value of the attracpotential and persistence length. We can collapse the vaof Dcoll on a single curve as a function of a proper combintion of the two dimensionless ratios,L/Lp and L0 /Lp ,(L/Lp)1/3(Lp /L0). Figure 6 shows the mastercurve for thdecay rate.
Table II reports the statistics of chains conformationsdifferent values ofL/Lp at two different times, i.e.,t1
6-6
-fote
eeitofu
taadnleo
teeth
ithlts.
rohde
ofto
ce.g
t, aslyedthehats ofweed,andA
s inost
ingre,herac-heterbleuch
ter-the
un-an-
u-
l-
COLLAPSE OF A SEMIFLEXIBLE POLYMER IN POOR . . . PHYSICAL REVIEW E69, 021916 ~2004!
5100t rot,Lpand t25200t rot,Lp
. At t1, the single head racquet is a conformation much more likely than the torusall the explored values ofL/Lp . This confirms that the direcformation of a torus is an unlikely event, as discussed in SIII A, and most of the chains collapse through a sequencfolding events. Att2, the percentage of molecules in th1-head racquet configuration does not vary substantially wrespect tot1, because while uncollapsed molecules foldform new 1-head racquets, some of these configurationsther collapse in higher order metastable racquets and stori. The number of chains collapsed in tori or multiple heracquets monotonically increases within the observed spatime. For shorter chains the likelihood of tori and multipracquets is similar, while the metastable states seem tinitially more favorable for longer chains (L/Lp58 andL/Lp510). Overall, the results show that the intermediaconformations in the collapse phenomenon are extremlong-lived: only 8.4% of the simulated chain has reachedequilibrium, torus state after 200t rot,Lp
.
D. Comparison of different computational models
We have performed simulations of collapsing chains wdifferent computational models in order to verify the resuobtained and the appropriateness of the model chosenparticular, we have performed simulations with a bead-model, where the chain connectors are rigid rods rather ta stiff elastic spring. As noted in Sec. II, the bead-rod mo
FIG. 5. Normalized collapse time distribution for all the simlation conditions.
02191
r
c.of
th
r-ble
of
be
slye
Indanl
is a constrained model which ensures local inextensibilitythe chain, but is computationally more expensive, duemore stringent time step requirements for convergenWhile local inextensibility is fundamental when analyzinthe polymer contribution to the stress tensor@22#, its effecton the pathways and dynamics of collapse is less relevanconfirmed by our results. In fact, we observe qualitativesimilar pathways to collapse in the simulations performwith the bead-rod model. The preferred pathway includesformation of metastable racquet shapes, similarly to wwas observed in the bead-spring simulations. The kineticthe collapse are slower for the locally inextensible model;believe that this is due to how the constraints are imposi.e., by projecting the random forces onto the constraints,therefore not allowing any fluctuations along the chain.more detailed comparison between the collapse kineticthe two models is prevented by the high computational cof simulating the constrained system.
We also perform simulations with a phantom bead-sprchain, i.e., a chain with no repulsive hard core. Once mowe find that the dynamics of collapse are qualitatively tsame, although the actual shapes of the tori and of thequets formed by the phantom chain are slightly different. Tabsence of hard core repulsion allows in fact for tighstructures, as shown in Fig. 7. In particular, in the metastamultiple racquets the neck for phantom chains appears mnarrower than for the chains where excluded volume inactions are taken into account, while it is noticeable thatshape and the curvature of the racquet head appearschanged. The traces in time of the energy level for the ph
FIG. 6. ~Color online! Masterplot of the decay rate of uncolapsed molecules.
7
TABLE II. Statistics of chain conformations for differentL/Lp at t15100t rot,Lpand t25200t rot,Lp
MONTESI, PASQUALI, AND MACKINTOSH PHYSICAL REVIEW E69, 021916 ~2004!
tom chain show the features already discussed in Sec. IIwell-defined energy states corresponding to metastable stures, with rapid transitions between them~see Fig. 7, inset!.However, the average life of metastable states appear sligshorter than in the complete model: intuitively, this is duethe wider range of motion of the phantom chain, which pmits folding paths otherwise prevented by the excluded vume interactions.
E. Effect of shear flow
We discuss here the results of some preliminary workthe effects of steady shear flow on the collapse dynamicsemiflexible chains. In particular we monitor the decay rDcoll at different flow strengths, and we compare the paways of collapse with those at equilibrium. Flexible chaiexpand in shear flow, i.e., the average end to end distaincreases while the molecules tend to orient with the flothe effect of the flow is therefore counteracting that of tattractive forces, which tend to form compact globuli. Stsemiflexible chains, withLp>L shrink in shear flow, due toa buckling instability@27#; such rigid chains, however, woulnot collapse to form tori and multiple racquets in bad svent, as the open, rodlike shape is stable at equilibrium@24#.In this work, we consider the collapse of less stiff semifleible molecules: is not cleara priori what to expect for thecollapse dynamics under shear flow.
We use a freely draining model for the chain under sheWhile hydrodynamic interactions are crucial in describithe dynamics of flexible chains in shear flow@28#, their ef-fect is far less relevant in a semiflexible chain. For flexib
FIG. 7. ~Color online! ~a! Metastable racquet-headed state fophantom chain and~b! corresponding correlation matrix. Inset i~a!: energy traces in time for collapsing phantom chains.
02191
:c-
tly
-l-
nofe-
ce;
-
-
r.
chains in fact, different segments of the chain contribute vdifferently to the viscous drag: the segments in the corethe polymer coil—or globule in the poor solvenconditions—are partially shielded from the flow, thus reduing the overall viscous drag. For a stiffer chain, however, teffect is negligible, because locally the molecule is straiand the shielding effect from the other part of the chainhighly reduced. The appropriate approach therefore is usslender body hydrodynamics; in this context, the friction tesor can be written asz5z iuu1z(I2uu). For a slender rod,z i50.5z. Because the anisotropy of the friction coefficieaffects the dynamics only marginally~although it affects theviscous component of the stress tensor!, in the present simu-lations we setz i5z ~isotropic drag! for simplicity.
We show here that strong enough shear flow consideraspeeds up the collapse kinetics because it increases thelihood that the two ends of the chain meet. Figure 8 shothe trend of the decay rate as a function of the Weissenbnumber Wi, where Wi5gt rot,Lp
, for three different values of
L/Lp . The Weissenberg number describes the ratio betwthe polymer characteristic relaxation time and the time scof the flow. Consistently with what was done in the othsections of this paper, we use the rotational diffusion timea rod of lengthLp as the polymer characteristic time. Foshear flow, the characteristic time scale is given by 1/g @29#.When the Weissenberg number Wi is smaller than onepolymer can relax—in this case rotate over a sectionLp—inthe time scale of the flow, while for Wi.1 the flow is fastwith respect to the polymer relaxation. At Wi50.1, the de-cay rate does not change much with respect to the equrium value, and there is not a clear trend among the differL/Lp . At Wi51, the increase in decay rate becomes instsignificant, as the average time to first contact halves;Wi510, the effect of shear flow is completely dominatinand the decay rate under the three different conditionsalmost equivalent. Under shear flow, the semiflexible mecules undergo a sequence of compressions and extenwhile tumbling in the plane of shear; in the compressiregion, the molecules tend to shrink, therefore raisingprobability that two sections far apart along the chain ba
FIG. 8. ~Color online! Plot of the decay rateDcoll vs Wi for
L/Lp53,5,8, u051kT, andRattr510/3.
6-8
thin
erfodheiguthsf t
inuroro
tet aina
oaywo
argiat
po
oftteretersthea
alar-e ofe-
a-ta-not
thelar,se
e is
ysngsed
ell-ibittestringdalultsorkiththees
they-in-
ol-
ttesthehisun-o.nceta-pu-
m
e
COLLAPSE OF A SEMIFLEXIBLE POLYMER IN POOR . . . PHYSICAL REVIEW E69, 021916 ~2004!
bone come in contact. The attractive forces then trapmolecule in the partially collapsed configuration, preventit from straightening back in the extension region.
However, the continuous tumbling of the chain and altnation of extension and compression does not allow themation of the higher order, compact structures observeequilibrium. In particular, the shear flow seems to inhibit tformation of the torus, which was never observed in the hshear flow simulations. Figure 9 shows two typical configrations observed in strong shear flow; while they showbasic features of the metastable racquets, the structureless compact and do not reach the same level of energy oequilibrium conditions.
IV. DISCUSSION
The dynamics of collapse of semiflexible moleculesbad solvent has been recently investigated in the literatDirect, visual experimental evidence of the formation of tand high order racquets has been obtained studying biopmers with sufficiently large persistence lengths@13,14#. The-oretical work and Monte Carlo simulations@11,12,24# haveconfirmed that the torus is the stable conformation, but inmediate, metastable racquet-shaped conformations exisare long lived. Previous Brownian Dynamics simulations2D @15# have firstly suggested the possible general pathwto the dynamics of collapse.
The present work provides new insights on the collapsesemiflexible chains, confirming the suggested pathwthrough 3D simulations, and showing the effect of the trelevant dimensionless ratiosL/Lp andL0 /Lp on the kineticsof collapse. The time series of the total energy of each chshow clearly the fast transition between well-defined enelevels, corresponding to different metastable intermedshapes, and the final, stable torus. We also showed thatpossible pathways towards the stable conformation are
FIG. 9. ~Color online! Collapsed configurations of semiflexiblchains in shear flow at Wi510; ~a! L/Lp58, and~b! L/Lp55.
er
02191
eg
-r-at
h-earehe
e.ily-
r-nd
ys
fs
inytewos-
sible: direct formation of a torus and successive foldingthe chain in progressively higher order racquets. The laone appears to be more probable for the range of paraminvestigated. The decay rate of uncollapsed states for alldifferent conditions of the simulation can be plotted asfunction of a proper combination ofL/Lp andL0 /Lp . Fur-ther work is needed to confirm the validity of this empiricscaling and verify it in a wider range of parameters. In pticular, this scaling does not provide a proper dependencthe decay rate from the persistence length alone, becausLpis never varied (Lp520 connector lengths for all the simulations!.
We also performed simulations with different computtional models, therefore confirming that the observed mestable conformations and pathways are general and dodepend on the model details. However, the kinetics, i.e.,actual value of the decay rate, can be different: in particufor a constrained bead-rod model the kinetics of collapappear slower, possibly due to the way the Brownian noisprojected onto the constraints.
We show preliminary results on the kinetics and pathwaof collapse of semiflexible chains in shear flow. Stroenough shear flow increases the decay rate of uncollapstates, but prevents the formation of the compact, wdefined configuration observed at rest, and seems to inhthe formation of the torus. Based on this evidence, the faskinetics of condensation should be achievable by sheathe molecules initially—to promote the initial collapse—anthen stopping the flow to let internal forces drive the fincollapse to the compact configurations. While these reshave been obtained in unbounded shear flow, future wwill consider the interaction of the collapsing molecules wa solid surface, which could be attractive or repulsive forchain. Also, a future possible direction of research includthe effect of local defects along the chain, which makemolecule locally more flexible or stiffer, as well as the dnamics of compaction when specific binding sites areduced in the chain, for example, by proteins on DNA mecules.
ACKNOWLEDGMENTS
The authors wish to thank Bernhard Schnurr and F. Gifor useful discussions, and R. Adam Horch for preparingvisualizations of the collapsing semiflexible molecules. Twork was partially supported by the National Science Fodation under Grant No. PHY99-07949, through award NCTS-CAREER-0134389, and through the Nanoscale Scieand Engineering Initiative award EEC-0118007. Computional resources were provided by the Rice Parallel Comtational Engineering Cluster~NSF-MRI-0116289! and by theSARA Computing and Networking Services in Amsterdathrough the Vrije Universiteit.
@1# P. G. de Gennes,Scaling Concepts in Polymer Physics, 1st ed.~Cornell University Press, Ithaca, 1979!.
@2# R. A. Bird, C. F. Curtiss, R. C. Armstrong, and O. HassagDynamics of Polymeric Liquids, 2nd ed.~Wiley, New York,
,
1987!, Vol. 2.@3# M. Doi and S. F. Edwards,The Theory of Polymer Dynamics,
1st ed.~Oxford University Press, Oxford, 1986!.@4# P.G. de Gennes, J. Phys.~France! Lett. 46, 639 ~1985!.
6-9
ica
m
, C
S
s
c
d.
iesnd. Ati-be
becedi.e.,
ted
MONTESI, PASQUALI, AND MACKINTOSH PHYSICAL REVIEW E69, 021916 ~2004!
@5# B. Ostrovsky and Y. Baryam, Europhys. Lett.25, 409 ~1994!.@6# K.A. Dawson, E.G. Timoshenko, and Y.A. Kuznetsov, Phys
A 236, 58 ~1997!.@7# A. Halperin and P. Goldbart, Phys. Rev. E61, 565 ~2000!.@8# B. Chu, Q.C. Ying, and A.Y. Grosberg, Macromolecules28,
391 ~2002!.@10# A.Y. Grosberg, Biofizika24, 32 ~1979!.@11# V.A. Bloomfield, Biopolymers44, 269 ~1997!.@12# H. Noguchi, S. Saito, S. Kidoaki, and K. Yoshikawa, Che
Phys. Lett.261, 527 ~1996!.@13# I. Baeza, P. Gariglio, L.M. Rangel, P. Chavez, L. Cervantes
Arguello, C. Wong, and C. Montan˜ez, Biochemistry26, 6387~1987!.
@15# F. Schnurr, F.C. MacKintosh, and D.R.M. Williams, EurophyLett. 51, 279 ~2000!.
@16# H. Noguchi and K. Yoshikawa, J. Chem. Phys.113, 854~2000!.
@17# T. Sakaue and K. Yoshikawa, J. Chem. Phys.117, 6323~2002!.@18# D.T. Hung, L.J. Marton, D.F. Deen, and R.H. Shafer, Scien
221, 368 ~1983!.@19# From slender body hydrodynamics,z54phs«, where «
[1/ln(L/r), r is the rod radius, andf («)'1 for slender rods.
02191
.
.
.
.
e
@20# A.Y. Grosberg, T.T. Nguyen, and B.I. Shklovskii, Rev. MoPhys.74, 329 ~2002!.
@21# P.S. Grassia and E.J. Hinch, J. Fluid Mech.308, 255 ~1996!.@22# M. Pasquali, V. Shankar, and D.C. Morse, Phys. Rev. E64,
020802~R! ~2001!.@23# M. Pasquali and D.C. Morse, J. Chem. Phys.116, 1834~2002!.@24# B. Schnurr, F. Gittes, and F.C. MacKintosh, Phys. Rev. E65,
061904~2002!.@25# See EPAPS Document No. E-PLEEE8-69-073402 for mov
showing two typical pathways, direct formation of a torus, asuccessive folding of the chain into higher order racquetsdirect link to this document may be found in the online arcle’s HTML reference section. The document may alsoreached via the EPAPS homepage~http://www.aip.org/pubservs/epaps.html! or from ftp.aip.org in the directory/epaps/. See the EPAPS homepage for more information.
@26# N. A. C. Cressie,Statistics for Spatial Data, 1st ed.~Wiley,New York, 1991!
@27# A. Montesi, C. Wiggins, and M. Pasquali~in preparation!.@28# D. Petera and M. Muthukumar, J. Chem. Phys.111, 7614
~1999!.@29# The behavior of semiflexible polymers in shear flow can
characterized also through the ratio of viscous forces induby the flow and elastic forces due to the bending stiffness,
the elasticity number El5gzL4/(a4kBTLp), wherea is theeigenvalue of the dominant buckling eigenmode. El is relato Wi through the ratioL/Lp .
