Date Published: September 7, 2012
Publisher: BioMed Central
Author(s): Iddo Aviram, Ilia Veltman, Alexander Churkin, Danny Barash.
Methods for simulating the kinetic folding of RNAs by numerically solving the chemical master equation have been developed since the late 90’s, notably the programs Kinfold and Treekin with Barriers that are available in the Vienna RNA package. Our goal is to formulate extensions to the algorithms used, starting from the Gillespie algorithm, that will allow numerical simulations of mid-size (~ 60–150 nt) RNA kinetics in some practical cases where numerous distributions of folding times are desired. These extensions can contribute to analyses and predictions of RNA folding in biologically significant problems.
By describing in a particular way the reduction of numerical simulations of RNA folding kinetics into the Gillespie stochastic simulation algorithm for chemical reactions, it is possible to formulate extensions to the basic algorithm that will exploit memoization and parallelism for efficient computations. These can be used to advance forward from the small examples demonstrated to larger examples of biological interest.
The implementation that is described and used for the Gillespie algorithm is freely available by contacting the authors, noting that the efficient procedures suggested may also be applicable along with Vienna’s Kinfold.
The RNA molecule, once considered as an intermediate step between DNA and proteins, has drawn much attention in recent years. Discoveries relating to its unique capabilities to prominently participate in gene regulation have motivated even more the concerted efforts to understand its folding and structural arrangement at several levels, both at the level of tertiary structure and that of secondary structure. The functional form of single stranded RNA molecules frequently requires a specific tertiary structure, but the scaffold for this structure is provided by secondary structural elements which are hydrogen bonds within the molecule. The four building blocks of RNAs are A,C,G,U and the base pairings among them form the secondary structure. This leads to several recognizable “domains” of secondary structure like hairpin loops, bulges and internal loops. Although the functional role of the RNA molecule in more detail is related to its three-dimensional structure, the RNA secondary structure is experimentally accessible and in many interesting cases may contain substantial important information to shed light on the relationship between structure and function. In general, RNA folding is thought to be hierarchical in nature [1,2], whereby a stable secondary structure forms first and subsequently there is a refinement to the tertiary fold. Thus, RNA conformational rearrangements that will be mentioned in the discussion can often be studied by examining their secondary structure, while keeping in mind the importance of tertiary structure.
The complete folding event is governed by the chemical master equation . In order to introduce the concept behind the reduction of the time-dependent RNA folding problem to that of stochastic chemical kinetics describing the time evolution of a well-stirred chemically reacting system, the Appendix follows closely references [22,23] in summarizing the formulation leading to Gillespie’s Stochastic Simulation Algorithm (SSA).
Here below, we demonstrate our SSA implementation on two “toy problem” examples. Our program that we have been developing is similar in style to Vienna’s Kinfold , but has been built in principle to have the ability of exploiting memoization and efficiency considerations as proposed above.
As a possible application of biological significance, the time-dependent approach discussed above is suggested for beneficial use in the problem of deleterious mutation prediction. To elaborate on this problem, a common way to detect deleterious mutations in the secondary structure of RNAs is to look for mutations that may cause a conformational rearrangement to occur. It was noted in  that there is some probability that even a single mutation can substantially alter the RNA secondary structure. Experimentally, this was observed in the spliced leader of Leptomonas collosoma, in RNA viruses [29,30], and in some other biological systems. Another very recent finding of biological importance is the existence of disease-associated Single Nucleotide Polymorphisms (SNPs) called “RiboSNitches” that have an RNA secondary structural consequence that results in a disease phenotype . Computationally, even before the added motivation as a consequence of the latter disease related finding, this gave rise to a procedure for detecting deleterious mutations using RNA folding predictions numerous times . Each time, relevant programs from an energy minimization package such as RNAfold from the Vienna RNA package [24,33] or Zuker’s mfold [34,35] can be used. In these packages, expanded energy rules  that were derived from an independent set of experiments are incorporated into the folding prediction algorithm. While the folding prediction problem described above is the most fundamental problem in RNA bioinformatics, the RNA mutation prediction problem is a sub-problem that uses the former multiple times, for various mutation combinations. Historically, initial works for the mutation prediction problem can be traced back to [37,38] and have been substantially revived in [32,39]. The first publicly available program for the RNA mutation prediction problem that takes into account only single-point mutation predictions was called RNAMute [40,41]. It uses the Vienna RNA package in its core. Subsequently, a web server dealing with similar issues was put forth called RDMAS . There are also some computationally challenging issues in the mutation prediction problem , mainly in the generalization to multiple-point mutations that can become computationally heavy if a ‘brute-force’ strategy of calculating all possible mutations is used without devising any unique approach. There have been various suggestions on how to reduce the number of mutations simulated or make the computations more efficient, for example [44-46]. In general, neither the original RNAMute  nor RDMAS  can handle multiple-point mutations. Consequently, RNAMute  was extended to MultiRNAMute  and based on the approach of , the web servers RNAmutants  and later corRna  were developed. A web server for MultiRNAMute was worked out in . There is, however, one common feature that should be taken into account when considering all of the programs dealing with RNA deleterious mutation predictions. Because several single point mutations inserted to the wildtype sequence can bring about to the same secondary structure, oftentimes there is a degeneracy in the mutations that is needed to be addressed. Any mutation prediction method for the purpose of conformational rearrangement in the secondary structure should therefore aim to report in each step (i.e., one-mutation, two-mutations, etc.) several mutation possibilities, not only a single one. If the method only reports a single possible mutation in each step, it probably means that a computational efficiency consideration was taken that may neglect many good candidate mutations that are conformationally rearranging just as well and will lead to the same secondary structure. Therefore, in order to fundamentally solve the degeneracy of mutations leading to the same structure, we suggest to perform for each one a time-dependent calculation and check how smooth the folding occurs in time, to discriminate those mutation candidates that get stuck in a local optimum for a while during the folding in time. This is quite an intensive calculation for sequences that are beyond “toy problems”, leading to a computational challenge from the numerical standpoint. It is also of considerable interest to check whether there is a correlation between the kinetic calculation and the static information obtained by performing energy minimization without taking into account what happens during the folding event. In order not to lose reliability, we suggest to consider all single point mutation combinations, and decide which one is the most likely to occur without getting trapped in a local minimum by using a time-dependent approach.
The significance of the initial development of efficient procedures described here can be divided into several items. First, the time-dependent folding is what takes place within the RNA molecule, and the static view of RNA structure only at the beginning and end may not be sufficient or complete in many cases. Experimental approaches to measure folding kinetics in detail, such as temperature jump experiments or single molecular methods , can be employed to check the computational model and predictions, in turn, can be pivotal to RNA rational design. Developing efficient numerical methods for the time-dependent folding simulation is therefore, by itself, an important goal. Here, we embarked on the stochastic approach, noting that if at all possible to achieve with practical computation time then one should definitely consider deterministic approaches  for the simulations of biologically relevant examples. Another direction for reducing computational cost is by an efficient exploration of discrete energy landscapes, which was developed in a recent work  by introducing a sampling method that allows for a fast yet accurate estimation of the transition probabilities between macrostates when coarse graining of the state space is used [4,16,50,52,53]. Second, a time-dependent approach to contribute in deleterious mutation prediction is suggested, which is still an open problem of considerable biological interest in a variety of RNA systems. For example, point mutations performed on an RNA virus such as HCV can alter virus replication  or lead to a dramatic reduction in translation initiation . Development of efficient time-dependent simulations can well assist from the predictive standpoint in such efforts.
The probability of a reaction Rkto occur in an infinitesimal time interval [t, t + dt) will be denoted by akXtdt known as the propensity function. Applying the law of total probability, one can obtain:
1. Evaluate akxtk=1M and asum=∑k=1Makxt.
The authors declare that they have no competing interests.
IA and IV and AC worked on the software design, carried out development and implementation, and participated in drafting the manuscript. DB conceived the study, coordinated the software design and drafted the manuscript. All authors read and approved the final manuscript.