# Glassy transition in a disordered model for the RNA secondary structure

## Abstract

We numerically study a disordered model for the RNA secondary structure and we find that it undergoes a phase transition, with a breaking of the replica symmetry in the low temperature region (like in spin glasses). Our results are based on the exact evaluation of the partition function.

The folded structure of biopolymers, like RNA and proteins, is crucial for understanding the biological functionality of these molecules  and its characterization still remains a challenging problem in statistical mechanics and theoretical biology . The folding problem usually consists in understanding if and how a particular biomolecule (maybe one selected by evolution and present now in nature) folds into its native conformation. In this Letter we are interested in the characterization of the most generic (i.e. random) RNA molecules. Even if real RNAs are not completely random, they present a very large variability in their sequences and no strong correlations in their bases. The interest in studying the limiting and somehow unphysical case of really random sequences arises in order to answer the following questions. Is the folding transition, that forces real biomolecules into their functional shapes, characteristic of those sequences selected by the evolution? Do random sequences show some phase transition too? We answer affirmatively to both questions, showing that the transition depends more on the geometrical constraints and on the interaction energies spread rather than on the specific sequence. However in the random case the transition is of a glassy type and the low-temperature phase is not dominated by a single native state. Our results may be very useful in order to understand better what could happen in a prebiotic world mainly made of random RNA sequences . Such a transition (partially found only in a very simplified model of proteins ) was suggested in previous studies of the RNA folding .

In this Letter we first study the thermodynamical properties of random RNAs, finding some hints for the existence of a glassy transition. The clear evidence for such a transition is shown in the last part of the paper and has been obtained thanks to the typical tools of disorder systems statistical mechanics: spin glass susceptibility and a related parameter (see Fig. 3). The connection with complex systems is well expected: the model has both disorder and frustration.

Generally speaking a classification among biopolymers includes a hierarchy of structures and in principle a complete description must include all these levels. RNA from this point of view is supposed to be simpler than DNA or proteins since its secondary structure seems to capture the essential features of the thermodynamics of the molecule. RNA molecules are linear chains consisting of a sequence of four different bases: adenine ($`A`$), cytosine ($`C`$), guanine ($`G`$) and uracil ($`U`$). The four bases are related by complementarity relations: $`CG`$ and $`AU`$ form stable base pairs with the formation of hydrogen bonds and are also known as Watson-Crick base pairs.

The secondary structure of RNA is the set of base pairs that occur in its three-dimensional structure. Let us define a sequence as $`\{r_1,r_2,\mathrm{},r_n\}`$, $`r_i`$ being the $`i^{th}`$ base and $`r_i\{A,C,G,U\}`$. A secondary structure on $``$ is now defined as a set $`𝒮`$ of $`(i,j)`$ pairs (with the convention that $`1ijn`$) according to the following rules:

a) $`ji4`$: this restriction permits flexibility of the chain in its three-dimensional arrangement.

b) Two different base pairs $`(i,j),(i^{},j^{})𝒮`$ if and only if (assuming with no loss of generality that $`i<i^{}`$):

$`i<j<i^{}<j^{}`$: the pair $`(i,j)`$ precedes $`(i^{},j^{})`$;

$`i<i^{}<j^{}<j`$: the pair $`(i,j)`$ includes $`(i^{},j^{})`$.

Condition b) avoids the formation of pseudo-knots on the structure and the resulting structure can be drawn on a plane. In real RNA structures it is known that pseudo-knots occur but are rare and they can be excluded as a first approximation .

The energy of a structure is simply defined as $`H[𝒮]=_{(i,j)𝒮}e(r_i,r_j)`$. Other phenomenological parameters (including stacking energies and loop penalties) could be considered in order to take into account the whole complexity of the energy function .

In our approach we assume a drastic approximation to the original model in order to improve its tractability both from numerical and analytical point of view. As a first step we consider sequences of only two symbols $`(A,B)`$, that appear with equal probabilities, and we assume that only two kind of base pairs occurs: $`AA`$ and $`BB`$ pairs with energy $`1`$ (in arbitrary units); $`AB`$ and $`BA`$ pairs with energy $`2`$. It is reasonable to assume that such a reduction of symbols will not affect the thermodynamical class of criticality of the model (this claim is supported by numerical results we have obtained with a 4-letter code and Watson-Crick base pairs). We did not remove the constraint which forbids the links on short distances, but we simplify it to: $`ji2`$. We think that this topological constraint must be kept in order to not drastically change the entropy of the model and then its thermodynamical behavior. In this model disorder (encoded in the sequence $``$) and frustration (induced by the planarity condition on $`𝒮`$) are clearly distinct. We hope this could make the model analytically more manageable.

The model can be formally considered as unidimensional and with long range interactions: the disorder giving rise to different interactions strengths, all with the same sign (here the disorder does not induce frustration), while the planarity condition making the long-distance links unlikely. We have numerically estimated that the probability of having a link between two bases a distance $`r`$ apart goes down roughly like $`r^{3/2}`$.

The planar structure of the configurations and the simple energy function chosen allow to write down  a recursion equation for the partition function of the subsequence contained inside the base interval $`(i,j)`$:

$$Z_{i,j}=Z_{i+1,j}+\underset{k=i+1}{\overset{j}{}}Z_{i+1,k1}e^{\beta e(i,k)}Z_{k+1,j},$$
(1)   
with $`Z_{i,i}=Z_{i,i1}=1i`$. Such a recursion relation is particularly effective since the time needed for the computation of $`Z_{1,L}`$ scales as $`O(L^3)`$. With a slight modification of the algorithm it is also possible to include similar recursions for the internal energy $`U=H[𝒮]`$ and its second moment $`U^{(2)}=H^2[𝒮]`$, where $``$ is the usual average over the Gibbs-Boltzmann distribution. At this level all the observables actually depend on the sequence over which they have been calculated and, if we want to gain information on the class of universality of the model, we have to average them over all the random realizations of the sequence.

In Fig. 1 we show the specific heat (averaged over the disorder) for sizes ranging from $`L=128`$ to $`L=1024`$. We note a very slow increasing of the peak height with the size, which seems not to diverge. There is no hint for a finite jump in $`C(T)`$. This could be compared with the result by Bundschuh and Hwa  who found a finite jump in the specific heat (note however that their model has an unique ground state, which dominates the frozen phase). It is important to point out that in the temperature region $`T0.150.2`$ the curves slightly cross themselves and as a consequence the decrease of $`C(T)`$ becomes steeper for larger sizes. One of the main effects of the disorder is that the location on the temperature axis of the critical region becomes sample-dependent. A measure of the critical region width can be achieved from the sample to sample fluctuation of the temperature where the specific heat has a peak ($`\mathrm{\Delta }T_p`$). We find that $`\mathrm{\Delta }T_pL^\omega `$, with $`\omega =0.26`$. If we assume that these fluctuations are induced by the presence of a nearby transition, we obtain a value $`\nu =\omega ^1=3.9(1)`$.

Since the model is unidimensional, $`\alpha =2d\nu =1.9(1)`$ and then the second derivative of the specific heat with respect to the temperature should display a very slow divergence or a finite jump. In fact, in the lower inset of Fig. 1, can be seen that the argument is fully supported by the data, which show the typical finite size behavior of a discontinuity. The clear crossing point of the data around $`T0.2`$ is supposed to be a signature for non-analyticities in thermodynamical potential. We note that such a point is located well below the peak temperature. This is a common feature in many disordered systems (e.g. spin glasses). Near this temperature also the entropy of the model has a crossing point, which signals a rapid shrink of the available phase space.

Moreover the model has a finite zero-temperature entropy (see upper inset in Fig. 1). The zero-temperature results have been obtain via an exact enumeration of all the ground states structures (GSS) for any given sequence. The number of GSS (i.e. the degeneracy) strongly depends on the sequence: for example, studying thousands of different sequences with length $`L=256`$, we have found sequences with degeneracies ranging from 1 to $`𝒪(10^7)`$. In the upper inset of Fig. 1 we show the zero-temperature entropy defined as $`S(T=0)=\mathrm{log}(𝒩)/L`$, where $`𝒩`$ is the GSS degeneracy and $`L`$ is the sequence length, as a function of $`L`$. The line is the power law extrapolation, which tends to $`S(T=0)=0.0255(8)`$ .

Since the model turns out to be highly degenerate in the low-temperature phase, the natural question is how these GSS are organized. It is quite obvious that a very different physical behavior may appear in a model whose GSS are all very similar (like an ordered or “ferromagnetic” behavior) compared to a model whose GSS are sparse over the whole configurational space. A more quantitative analysis can be achieved introducing the notion of distance between structures and a classification based on these distances. To quantify the relative distance between two structures, we have used the overlap, which is defined as

$$q[𝒮,𝒮^{}]=\frac{1}{L}\underset{i<j}{}l_{ij}^{(𝒮)}l_{ij}^{(𝒮^{})},$$
(2)   
where the variable $`l_{ij}^{(𝒮)}`$ ($`l_{ij}^{(𝒮^{})}`$) takes value 1 if sites $`i`$ and $`j`$ are connected in the $`𝒮`$ ($`𝒮^{}`$) sequence and 0 otherwise. By definition the overlap takes values in the interval $`[0,1]`$. For any given disorder realization (i.e. sequence) $``$ we can define the zero-temperature probability distribution function (pdf) of the overlaps as

$$P_{}(q)=\underset{𝒮,𝒮^{}\mathrm{\Gamma }_{}}{}\delta (qq[𝒮,𝒮^{}]),$$
(3)   
being $`\mathrm{\Gamma }_{}`$ the GSS set. This definition can be easily generalized to every temperature summing over all the structures and weighting each term with the Gibbs-Boltzmann factor of $`𝒮`$ and $`𝒮^{}`$. The usual classification of disordered systems  is based upon the average pdf of the overlaps, the so-called $`P(q)[P_{}(q)]`$, the average being taken over the disorder distribution function.

We have calculated the $`P(q)`$ at different temperatures, $`T[0,0.4]`$. While at $`T=0`$ we summed over the whole sets $`\mathrm{\Gamma }_{}`$, at finite temperatures we performed a Monte Carlo sampling of the structures in the spirit of Higgs .

In Fig. 2 the averaged $`P(q)`$ are shown. The first striking evidence is that, decreasing the temperature, the shape of the $`P(q)`$ changes abruptly from a narrow peak in the low-$`q`$ region to a broader one which extends over almost the whole allowed support. In the insets we present the size dependence of $`P(q)`$ for the highest and lowest temperature considered. For the $`T=0.4`$ case, we are highly confident that the thermodynamical limit would be a delta function (the width of the distribution goes to zero as $`\mathrm{\Delta }qL^{1/2}`$). For the $`T=0`$ case, the asymptotic shape is much more difficult to be extrapolated, since the width of the $`P(q)`$ scales with a small power of $`L`$ (as in ) and, eventually, we can not exclude that it goes to a finite value, implying a breaking of the replica symmetry.

While in Fig. 2 the averaged $`P(q)`$ gives us information about the typical pdf of the overlaps, we can get some hints about the origin of the $`P(q)`$ broadness in the low-temperature phase analyzing directly the $`P_{}(q)`$ for each sequence. If all the GSS of a given sample are very similar its $`P_{}(q)`$ will be non-zero only in a narrow $`q`$-range not too far from the upper bound $`q=1`$. On the other hand, if the GSS are very heterogeneous their mutual overlaps will cover a large $`q`$-range.

The great majority of the sequences shows a very broad $`P_{}(q)`$, signaling a strong heterogeneity in the GSS. Moreover the shape of the pdf completely changes from sequence to sequence (this property is called non-self-averageness in spin glass jargon ). Nevertheless some patterns can be easily recognized: while single peak shapes are mostly associated with low-degeneracy sequences, highly structured ones seem to be not correlated to their degeneracy and they are responsible for the $`P(q)`$ broadness. Among the latter the double-peak shape dominates, especially for the sequences with higher entropy: the higher $`q`$ peak gives information about the typical distance between two structures in the same state , while the lower $`q`$ one can be associated with the rising of a backbone , that is the set of persistent links common to all the GSS (already found in ). The position of this second peak strongly fluctuates from sample to sample, giving rise to the long tail in the $`P(q)`$, like spin glass models in external field.

In order to understand whether a true transition happens in this model, we have measured the order parameter introduced in 

$$A=\frac{[\chi _{}^2][\chi _{}]^2}{[\chi _{}]^2},$$
(4)   
where $`\chi _{}=L(\mathrm{\Delta }q_{})^2`$, being $`\mathrm{\Delta }q_{}`$ the width of $`P_{}(q)`$. The $`A`$ parameter measures how much the $`P_{}(q)`$ changes from sample to sample. The crossing point of different curves in Fig. 3 signals the existence of a low-temperature spin glass phase, where the $`P_{}(q)`$ become non-self-averaging (analogous results has been obtained with the 4-letter model). In this phase the RNA is “folded”, that is the number of links is nearly the maximum allowed.

The critical temperature of this transition seems to be located between $`T=0.1`$ and $`T=0.15`$. We have determined the best estimates for $`T_c`$ and for the critical exponent $`\eta `$ requiring the best collapse for the susceptibility $`\chi [\chi _{}]`$ data, scaled assuming the usual finite-size formula $`\chi =L^{2\eta }f(L^{1/\nu }(TT_c))`$ (see inset of Fig. 3). We obtain the values $`T_c0.13`$ and $`\eta 1.41`$. We stress that we also tried to collapse the data fixing $`T_c=0`$, but the result was very poor.

The critical temperature seems to be below the one we found by the study of thermodynamical quantities. However given the high value of $`\nu `$, the critical region should shrink as $`L^{1/\nu }`$ and then all the region around $`T=0.10.3`$ is critical as suggested by the wide separation of the two peaks in $`_T^2C`$ (lower inset of Fig. 1).

We have presented strong numerical evidences for a phase transition in a random model for the RNA secondary structure. It is very important to stress that the thermodynamical limit is not so interesting for biological RNAs, which are at most thousands bases long. As a consequence our sizes are in principle directly comparable with a large number of biological molecules. Our findings about the broadness of the $`P(q)`$ could suggest for the existence of zero-energy fluctuations of the order of the volume, which is a well known behavior in spin glass and disordered systems. In , for example, it has been found that the matching problem (which is disordered and frustrated) has low-energy excitations of order $`\sqrt{L}`$. These excitations becomes irrelevant in the thermodynamical limit, but they are a key ingredient in order to correctly describe finite systems. In the low temperature region of our model (from $`T=0.13`$ down to $`T=0`$) $`\chi L^{0.6}`$ and it has strong fluctuations from sample to sample. This situation can be described by an effective breaking of the replica symmetry with a strength which goes to zero as $`L^{0.4}`$ according to the slow shrinking of $`P(q)`$ at $`T=0`$. Incidentally we note that in all this temperature region the critical exponent of $`\chi `$ is the same, as suggested from the scaling plot in the inset of Fig. 3. Moreover we have also measured the $`G`$ cumulants defined in  and we have verified that it goes to the value $`\frac{1}{3}`$ as the temperature goes to zero coherently with a replica symmetry breaking scenario.

In conclusions, we have found a glassy transition in a simplified random model for the RNA secondary structures. This transition corresponds to the breaking of the configurational space in many disconnected regions (ergodicity breaking). In terms of random RNA folding, this means that below the critical temperature almost every sequence folds (all the low-energy structures are very compact), but very often not in a single structure. The ergodicity breaking is of primary importance also for the folding dynamics, that may become very slow (glassy).

We have checked that the transition disappears as soon as we remove the constraint of not having links on short distances (maybe this is a pathology of the 2-letters code) or as soon as we set all the interaction energies to the same value. These facts suggest us that the glassy transition is mainly due to the freezing of some strong links, which then force the rest of the interactions, aided by the geometrical constraints. A cooperation phenomenon between interaction energies heterogeneity and geometrical constraints has been already observed in DNA models .

We warmly thank R. Zecchina for many interesting discussions and for a careful reading of the manuscript.