Model for Folding and Aggregation in RNA Secondary Structures

0.23 MB·English

by  Vishwesha Guttal

A Model for Folding and Aggregation in RNA Secondary Structures
Vishwesha Guttal and Ralf Bundschuh
191 W Woodruff Ave, Department of Physics, The Ohio State University, Columbus, Ohio 43210-1117
We study the statistical mechanics of RNA secondary structures designed to have an attraction
between two different types of structures as a model system for heteropolymer aggregation. The
competition between the branching entropy of the secondary structure and the energy gained by
pairing drives theRNA toundergo a ‘temperature independent’ second order phase transition from
6 a molten to an aggregated phase. The aggregated phase thus obtained has a macroscopically large
0 numberofcontactsbetweendifferentRNAs. Thepartitionfunctionscaling exponentforthisphase
0 is θ ≈ 1/2 and the crossover exponent of the phase transition is ν ≈ 5/3. The relevance of these
2
calculations to theaggregation of biological molecules is discussed.

n
a PACSnumbers: 87.15.Aa,87.15.Cc,64.60.Fr,87.15Nn
J

2
1 RNA secondary structures are an excellent model sys- aggregation is associated with various diseases such as
tem to study folding phenomena in heteropolymers. Un- Alzheimer’s, Mad Cow, and Parkinson’s [8, 9]. Thus,

] like in the protein folding problem where a large num- this phenomenon has been studied with the tools avail-
h
ber of different monomers needs to be taken into ac- able for the protein folding problem in various con-
p
count to understand folding [1], an RNA has just four texts [9, 10]. But also in the realm of RNA folding the
-
o bases A, U, C, and G. The interactions between these competition between individual and aggregated struc-
i
b bases are simpler than in the protein folding problem turesplaysanimportantrole,e.g.,inthegrowingfieldof
. due to the separable energy scales of the secondary riboswitches [11]. In these riboswitches, the aggregation
s
c and the tertiary structure. These features make RNA of two RNA molecules through base-pairing in competi-
i secondary structures a both analytically and numeri- tion with the base pairing of the individual molecules is
s
y cally amenable model for rigorously studying various used to regulate the expression of genes in dependence
h genericthermodynamicpropertiesofheteropolymerfold- on the concentrations of the RNAs involved. Even the
p ing [2, 3, 4, 5, 6, 7]. local structure of double stranded DNA in the repeat
[ Quite a lot is known about the folding thermodynam- regions of the genes involved in triplet repeat deseases

2 icsofsingleRNAmolecules. Atlowtemperatures,where (Huntingtons, fragile X, etc. [12]) is an example of an
v monomer specific binding energies and sequence hetero- aggregated structure (the double stranded DNA) com-
0 geneity are important, the resulting (frozen) phase is peting with the multitude of secondary structures the
4 glassy [3, 4]. At high temperatures, thermal fluctuations single strands of this DNA can form by themselves since
1
lead to a denatured phase, where the backbone is ran- their repeat units of, e.g., CAG and CTG in Hungting-
4
domly coiled (without any binding) like a self-avoiding ton’s desease, allow self-pairing as well. Here, we ap-
0
random walk. At intermediate temperatures, where an proach the phenomena associated with competition be-
5
0 effective attraction between short segments is impor- tween intra-molecular structure and aggregationby con-
/ tant, the molecules are expected to be in the so called sidering a toy model to study the phase transition of an
s
c molten phase [4, 5]. In this molten phase many differ- RNA secondary structure from the molten to an aggre-
i ent secondary structures all having comparable energies gated phase. While our model is literally applicable to
s
y (within O(kBT)) coexist. If the tendency of biological the abovementioned triplet repeatdesease genes,we see
h sequences to be designed to fold into a specific, func- itmorebroadlyasabasicmodelforstudyingthecompe-
p tional structure is taken into account, the native phase titionbetweenintra-molecularstructureandaggregation
v: emerges [5, 6]. Many important questions have been into which later aspects of the other scenarios discussed
i raised with regard to these phases, e.g. their stability, abovesuchasnativestatesandsimultaneousaggregation
X
characteristics, and the properties of the phase transi- of several molecules can be incorporated. In studying
r tions between them in the context of both protein and ourmodel,ourfocusisonthethermodynamicproperties
a
RNA folding [1, 2, 3, 4, 5, 6, 7]. In this Letter we of the system. Thus, we solve the model exactly in the
shall begin to understand another important aspect of thermodynamiclimitandcalculatethecriticalexponents
heteropolymer folding, namely the competiton between relevant to the phase transition.
the individual folding of the molecules and aggregation RNA is a biopolymer with four different monomers A,
of severalmolecules, using the RNA secondarystructure U,CandGinitssequence. TheWatson-CrickpairsA-U
formulation. andC-G are energeticallythe most favorablepairs while
In the context of protein folding, the competition be- G-U is marginallystable andthe other combinationsare
tween individual folding and misfolding associated with prohibited. By an RNA secondary structure, we mean
aggregation is a very important phenomenon. The fail- a collection of binding pairs (i,j) with 1 i < j N,
≤ ≤
ure of protein molecules to fold correctly and the asso- where N is the number of bases in the sequence. Any
ciated formation of alternative structures stabilized by two pairs (i1,j1) and (i2,j2) are either nested, i.e. i1 <
2
beofcomparablemagnitudeinarealisticRNAmolecule.
The inter RNA base pairing energy, or the bias, ǫ3 is
the parameter which can in principle be controlled by
sequence mutation. Note, that neglecting sequence het-
erogeneity in this kind of models was established as a
useful approximation at not too low temperatures in a
similar context [5] (see also the discussion at the end of
this letter).
FIG. 1: Abstract representations of RNA secondary struc- Denote the Boltzmann factors corresponding to the
tures (from [4]). (a) Helix representation (b) Non-crossing pairing energies by q1, q2 and q3 respectively. We show
Arch diagram. Here, the solid line corresponds to the back- that this simple model predicts a molten to an aggre-
bone of the RNA.The dashed arches correspond to the base gatedphase transition,as wetune the parameterq3. We
pairs. The absence of pseudo-knots implies that the arches do so by exploiting the recursive relation [14, 15]
never cross. (c) Mountain Representation. Here, as we go
along the backbone of the RNA from base 1 to N (repre- j−1
sented bythebase line), wego onestep upfor thebeginning Zij =Zi,j−1+ Zi,k−1e−ǫjk/TZk+1,j−1 (2)
of a pair, one step down for the closing of a pair and a hori- X
k=i
zontalstepfornopairing. Suchamountainnevercrossesthe
baseline and always returnsto thebaseline at theend. forthepartitionfunctionZ forasequenceofbasesfrom
ij
i to j, which can be evaluated in O(N3) time starting
i2 < j2 < j1 or independent i.e. i1 < j1 < i2 < j2. froTmotkheeepintihtiealacnoanlydtiitcioanlscaZlci,uil=atiZoni,is−s1im=p1le., we assume
Theaboverestrictionmeanswearenotallowingpseudo-
each RNA to be of equal length, containing N 1 bases
knots, which are generally energetically not as favorable −
[16]. WenowconsiderthejointfoldingofthesetwoRNAs
[13]. Sucha secondarystructurecanbe representedby a
helix diagram, non-crossing arch diagram or a mountain anddenote its partitionfunction by Zd(N;q1,q2,q3). As
explained before, the free energy of pairing for the bases
representation as shown in Fig. 1.
belonging to a given RNA has contributions from the
Letthefreeenergyassociatedwiththepairingofbases
energy gain due to the pairing and the entropy loss as-
i and j in an RNA be ǫ . This free energy has contri-
ij sociated with the loop formation. This holds true even
butions from the gain in the energy due to binding and
forpairing acrossthe bases belongingto different RNAs.
the associated configurational entropy loss. In addition
Butwhenthefirstpairingbetweenthebasesbelongingto
to these, in principle there are also entropic and/or en-
differentRNAsoccur,thereisanadditionalentropicloss
ergetic effects due to loopformation,stacking,etc. Even
due to the breakdown of translational invariance sym-
thoughtheaccurateparametersasdeterminedbytheex-
metry. Thereafter,only the free energy ǫ3 plays a role in
periments are essential to calculate the exact secondary
theinterRNAbasepairing. Inthethermodynamiclimit,
structure, such microscopic details as well as the ex-
this additionalentropiclosshasnoeffect onthe phaseof
act values of the energies ǫ do not affect the asymp-
ij the system, but it is the energetics of pairing that drives
totic properties of the phases and the critical exponents.
the phase transition. Hence, we ignore this additional
Hence, we ignore them in our model calculations.
entropic term. This essentially reduces the problem to
To understand the phase transition from the molten
the folding of a single sequence with 2N 2 bases. The
to the aggregated phase, we first define the aggregated −
aggregated secondary structure can now be interpreted
phase as an ensemble of RNA secondary structures in
ashavingamacroscopicallylargenumberofcontactsbe-
which a macroscopically large number of contacts occur
tween the two halves of the concatenated RNA.
between two different RNAs. We consider a dual RNA
Let us first consider two special cases. Setting q =
biomolecule system consisting of two types of RNA in a
q1 =q2 =q3 correspondstothewellknownmoltenphase
solution. We refer to them as RNA-1 and RNA-2. In-
of the RNA secondary structure, whose partition func-
dividually, RNA-1 and RNA-2 are in the molten phase.
tion can be calculated exactly in the asymptotic form
Hpaoiwrienvgesr,bwethweenenthbeaysaesrefrtoomgetdhiffereriennat msooluletciounle,saalsroebpaosse- Zd(N;q,q,q) = Z0(2N;q) = A(q)(2N)−θzc(q)2N with
the characteristic scaling exponent θ = 3/2 [5]. This
sible. We study the phases of this dual RNA system, as
exponent is characteristic in the sense that it is insen-
the bias strength is varied.
sitive to various microscopic details of the RNA sec-
To do so, we assume a simple pairing energy model
ondary structure such as the cost of a hairpin loop,
with the free energy of pairing between bases i and j
weak sequence heterogeneity, etc. The other simple
defined as:
case is q3 = 0. This case describes two RNAs in the
ǫ1 if i,j RNA-1 molten phase which do not know of each other’s pres-
ǫi,j = ǫ2 if i,j ∈RNA-2 (1) ence. The partition function of such a dual RNA is then
 ǫ3 if i ∈RNA-1,j RNA-2 or vice-versa just the product of individual partition functions, i.e.
 ∈ ∈ Zd(N;q1,q2,0) ≡ Z0(N,q1)Z0(N,q2). Hence the scaling
Here,theintraRNAbasepairingenergiesǫ1andǫ2could exponent is θ =3.
3
q3 = 0 in the asymptotic limit, hence θ = 3. This claim
is verified by numerical calculations of the exact parti-
tion function for finite length and the calculation of an
asymptotic macroscopic quantity (the order parameter)
tobe definedbelow. Theresultingsimple phasediagram
is shown in the inset of Fig. 2.
Inordertoverifythatthephasetransitionindeedhap-
pensatq3c =√q1q2,wecalculatethe orderparameterof
the phase transition. Here, the order parameter Q is de-
finedasthefractionofsecondarypairingsinasecondary
structure,animportantstructuralpropertyoftheaggre-
gate. For arbitrary q3 the order parameter can be cal-
culated exactly from Q = limN→∞(dlnZd/dlnq3)/N.
−
The inset of Fig. 3 clearly shows Q = 0 for q3 q3c
≤
FIG. 2: (color online) The behavior of the partition function and continuously increasing with q3 thereafter saturat-
Zˆd(z;q1 = 4,q2 = 9,q3). For q3 = q3c = 6, we observe a ing to Q = 1 for q3/q3c 1. From this behavior of the
≫
square root behavior. For q3 >q3c, we see an inverse square order parameter we can conclude that the phase transi-
root behavior. The inset shows the resulting phase diagram. tion indeed occurs at q3c = √q1q2 and that the phase
transition is of second order. Physically, we can under-
stand the behavior of the order parameter by using the
We now want to understand the case of general q1, mountain representation of RNA (see Fig. 1c). Between
q2 and q3. To this end we calculate the partition func- anytwoconsecutivesecondarypairings,thecontribution
tion of the dual RNA as follows. Let the base pairings of primary pairs to the height of the mountain is zero.
within a given RNA be called primary and those across Hence,thetotalnumberofsecondarypairingsisequalto
differentRNAsbecalledsecondary. Anygivensecondary theheight h ofthemountainatitsmidpoint. Usingthe
structurethusobtainedhasaseriesofsecondarypairings randomwahlkianalogy[4,18],wefindthat h O(N1/2),
(i1,j1),...,(ik,jk)suchthat1≤i1 <...<ik ≤N−1and hence Q ∼ O(N−1/2). For q3 < q3c, the shecio∼ndary pair-
1 j1 < ... < jk N 1. Note that we have labeled ingsareevenlesslikely,andhenceinthethermodynamic
≤ ≤ −
the RNA-1 by i and the RNA-2 by j indices. The bub- limit Q = 0 for q3 q3c, consistent with what we have
≤
bles thus formedbetweenany twoconsecutivesecondary obtained by exact expression.
pairings are allowed to have only the primary pairings.
To further verify our claims about the phase for q3 <
If all the secondarystructure configurations are enumer-
q3c and to calculate the scaling exponents correspond-
ated according to the number of the inter-RNA (or the
ing to the second order phase transition, we iterated the
secondary) contacts k, then the total partition function
recursion relation (Eq.(1)) to calculate the exact parti-
of this dual RNA system, in the z-transformrepresenta-
tion function for RNA of finite length N. The results
tion can be written as:
of the numerical calculations are in complete agreement
∞ with the phase diagram of Fig. 2 (inset) when extrap-
Zˆd(z;q1,q2,q3)= q3kZˆ0(z;q1)k+1 Zˆ0(z;q2)k+1 (3) olated to the thermodynamic limit, thus verifying our
X ∗ claim [19]. Next we calculate the free energy per length
k=0
dz′ Zˆ0(z′;q1)Zˆ0(z/z′;q2) f(q1,q2,q3) = −lnZd(N)/N, taking into account the fi-
= (4) nite size effects. We assume the usual scaling function
I z′ 1−q3Zˆ0(z′;q1)Zˆ0(z/z′;q2) for the order parameter Q(N)=N−1/2g[(q3−q3c)N1/ν]
closetothe criticalpoint. Fig.3showsthe resultofscal-
whereZˆd(z;q1,q2,q3)andZˆ0(z;q)arethez-transformsof ing plot, with the best fit value for the crossover critical
Zd(N;q1,q2,q3) and Z0(N;q) respectively. The symbol exponent ν 5/3.
indicates the convolution in z-space defined as f g = ≈
∗ dz′f(z′)g(z/z′). Eq.(4)isobtainedbysumming u∗pthe This model has some similarities with G¯o-like model
z′ studiedbyBundschuhandHwa[5]whichshowsamolten-
Hgeometric series in Eq.(3). The convolution integration
native transition. The physics behind the phase transi-
canbe done numericallytoobtainthe singularitiesofZˆd tion in their model as well as our model is the same,
and hence, the asymptotic behavior of Zd(N;q1,q2,q3). i.e., the competition between the energetic gain of the
The resultsareshowninFig.2. Forq3 =q3c =√q1q2, secondarycontacts (or native contacts of G¯o-like model)
we find a square root singularity and hence θ = 3/2 and the branching entropy. But, contrary to the native
[17], the characteristic exponent of the molten phase. phase where the ground state is unique, the aggregated
For q3 > q3c, Zˆd has an inverse square root singular- phase has degenerate ground states. On the other hand,
ity, indicating a new phase. We interpret the new phase both these models can ‘melt’ from their (aggregated or
with the partition function scaling exponent θ 1/2 as native) ground state to any of the molten, glassy or de-
≈
the aggregated phase. We claim that for all q3 < q3c, natured phase, depending on the temperature and the
the dual RNA system is just the phase corresponding to strength of the bias. The differences in the behavior of
4
work [4, 5, 19]. However, at low enough temperatures
1
2 0.8 q1=4 and q3=9 where such a homogeneous approximation is no longer
valid, it should be interesting to consider the role of a
0.6
Q suitably defined bias in the glassy phase. Secondly, we
2 0.4
1/N 0.2 considered only two RNA molecules for aggregation. In
N) 0 the case of multiple RNAs participating in the folding,
Q(1 0 2 q4/q6 8 10 N=100 thegroundstatewoulddependonhowthedifferenttypes
3 3c N=400
N=200 of RNAs are alignedto fold. In fact, these groundstates
N=800 could be topologically different from the ground state of
N=1600
our two RNA model. Hence, the values of the critical
−050 −25 0 25 50 75 100 exponents for the transition might change, though the
(q - q )N1/ν qualitative physics of aggregation, such as the critical
3 3c inter-molecular base pairing energy at which the transi-
tion takes place would remain the same.
FIG. 3: (color online) Scaling plot for the order parameter.
Inset shows the order parameter of the phase transition. In Insummary,wehavepresentedasimplemodelforhet-
both the plots q1 =4 and q2 =9, henceq3c =6. eropolymer folding using the RNA secondary structure
formulation, which shows a second order phase transi-
tionfromanindependentlymoltentoanaggregatedphase.
these models arises from the fact that for the G¯o-like
The behavior at criticality turns out to be the molten
model the bias is site specific where as for the model
phase for the concatenated molecule. The transition is
we havepresented, the bias is towards a macroscopically
completelydrivenbytheenergeticsofpairingandistem-
large number of sites.
perature independent. Proteins are known to undergo a
We would like to emphasize two simplistic consider-
foldingtransitionfromthenativetoanaggregatedphase
ations of our model. Firstly, choosing a molten phase
instead of from a molten to an aggregated phase [9]. It
to start with, in which the sequence heterogeneity is
shouldbeinterestingtoseeifthis studycanbeextended
unimportant. It is important to note that the bases as
tounderstandthe thermodynamicsofsuchaphasetran-
we call them here are not necessarily single bases, but
sition. It should also be interesting to study the role of
short segments of the sequence (such as CAG) whose ef-
kinetics of RNA folding in this phase transition.
fective interaction with any other segment is the same
[5]. Even if we do have weak sequence heterogeneity, We gratefully acknowledge useful discussions with
we do not expect that such microscopic details alter the Tsunglin Liu. RB is supported by the National Science
thermodynamicresultspresentedherebasedonprevious Foundation through grant no. DMR-0404615.
[1] K.A.Dilletal.,ProteinSci.4,561(1995);J.N.Onuchic Prusiner (Cold Spring Harbor Laboratory Press, Cold
at al, Annu.Rev. Phys. Chem. 48, 545 (1997); T. Garel Spring Harbor NY,1999), p. 1.
et al., J. Phys. I 7, 27 (1997); E. I. Shakhnovich, Curr. [9] forareviewsee: C.M.Dobson,Nature,426,884(2003).
Opin.Struct.Biol. 729(1997). [10] A. Slepoy, R.R.P. Singh, F. Pazmandi, R.V. Kulkarni,
[2] P.G. Higgs, Q.Rev.BioPhys. 33, 199 (2000). D.L. Cox , Phys. Rev.Lett. 87, 058101 (2001).
[3] P.G.Higgs,Phys.Rev.Lett.76,704(1996);A.Pagnani, [11] G. Storz, S. Altuvia, and K. Wassarman, Annu. Rev.
G. Parisi and F. Ricci-Tersenghi, Phys. Rev. Lett. 84, Biochem. 74, 199 (2005).
2026(2000);A.K.Hartmann,Phys.Rev.Lett.86,1382 [12] M. Mitas, Nucleic AcidsRes. 25, 2245 (1997).
(2001); A. Pagnani, G. Parisi and F. Ricci-Tersenghi, [13] I. Tinoco Jr. and C. Bustamante, J. Mol. Biol. 293, 271
Phys.Rev.Lett.86,1383(2001);F.Krzakala,M.M´ezard (1999) and references therein.
andM.Mu¨ller, Europhys.Lett.57,752(2002); E.Mari- [14] J. S.McCaskill, Biopolymers 29, 1105 (1990).
nari, A. Pagnani and F. Ricci-Tersenghi, Phys. Rev. E. [15] M. S.Waterman,Advances in Mathematics, Supplemen-
65, 041919 (2002). tarystudies,editedbyG.-C.Rota(Academic,NewYork,
[4] R. Bundschuh and T. Hwa, Phys. Rev. E, 65, 031903 1978), pp.167-212.
(2002);R. Bundschuh and T. Hwa, Europhys. Lett. 59 [16] Theequallengthapproximationisnotnecessary.Ingen-
903 (2002). eralaslongastheRNAsareoftheorderofsamelength,
[5] R. Bundschuh and T. Hwa, Phys. Rev. Lett. 83, 1479 all thesubsequent results hold.
(1999). [17] IftheZ-transformofthepartitionfunctionshowsapower
[6] P.G. Higgs, J. Phys. I (France), 3, 43 (1993). law singularity at zc, say Zˆ(z) ∼ (z −zc)α, then the
[7] H.OrlandandZ.Lee,Nucl.Phys.B,620,456(2002);R. corresponding partition function scaling exponent θ =
Mukhopadhyay et al., Phys. Rev. E. 68, 041904 (2003); α + 1. See App. A of Ref T. Liu and R. Bundschuh,
M.Baiesietal.Phys.Rev.Let91,198102(2003);P.Leoni Phys. Rev.E, 69, 061912 (2004) for derivation.
and C. Vanderzande,Phys. Rev.E. 68, 051904 (2003). [18] See,eg.W.Feller, An Introduction to Probability and its
[8] S.B. Prusiner, in Prion Biology and Diseases, ed. S.B. Applications (Wiley, NewYork,1950).
5
[19] Similar numerical calculations for the repetitive se- agree with these results.
quences with monomer specific binding energies also