Computational Design of a Circular RNA with Prionlike Behavior

During the last decade, the field of synthetic biology has impressively illustrated that nucleic acids and in particular RNA molecules are reliable materials for the design and implementation of functional circuits as well as nano-scale devices and objects [18, 24, 1, 23]. The reasons for this success are grounded in the facts that for RNA (i) an experimentally measured energy model exists, (ii) regulation at the level of RNA molecules is faster than via the production of proteins, and (iii) design questions are more readily expressed in the discrete framework of binary base pairing than in continuous interactions between, say, the amino acids in proteins.

RNA molecules have been extensively engineered in the classical context of gene regulation. Successful designs include Boolean networks with miRNAs [27, 34, 47], synthetic RNA switches [22, 40, 17], and artificial ribozymes [43, 13]. In general, the RNA is optimized to switch conformations upon a hybridization interaction, in order to toggle between an on and an off state. Turberfield et al. [38] showed for DNA that such hybridization reactions are reversible by a mechanism called toehold exchange, which was later used to design reversible logic circuits [14] as well as transcriptional oscillators [25] and transcription regulators [37]. Recently, multiple toehold switches have been integrated in vivo to activate gene expression in response to endogenous RNAs [17]. Cascading of strand-displacement reactions allows for the construction of multi-component chemical reaction networks, exhibiting complex computational and information-processing abilities [50, 41].

DNA and RNA hybridization effects are essential features for algorithmic self-assembly [36, 48] in nanotechnology. RNA nanoparticles have been constructed from smaller self-assembling units [6, 5, 16]. Small RNA motifs, such as complementary kissing hairpins, can be used to assemble complex 1D, 2D, and 3D shapes (for a recent review see [15]). Cayrol et al. [5] have found natural self-assemblies of the small RNA DsrA in E. coli, suggesting concentration-dependent RNA regulatory mechanisms via self-assembly.

The inherent complexity of nucleic-acid-based self-assembling systems makes it necessary to optimize hand-crafted designs computationally. Several energy-directed computational methods have been devised for the rational design of nucleic acid molecules that fold into single [20] or multiple [11, 21] predefined conformations or that form ensembles of interacting nucleic acid strands [49, 44]. While these methods carefully model equilibrium properties of the designed RNAs, they provide rudimentary or no support for designing kinetic properties, such as refolding times, which are much more expensive to compute and thus remain a challenge for computational design methods.

The protein-only hypothesis for the scrapie agent (for a review see [2]) proposes that a prion protein, with an altered (infectious) β-sheet-rich conformation, starts an autocatalytic cascade that uses the normally folded prion proteins as a substrate, converting them to the infectious form. This altered conformation then either self-assembles into fibers, which is the usual phenotype upon prion infection, or catalyzes the refolding of the remaining normally folded prions. A high activation energy between the normal and the infectious conformation prevents spontaneous conversion at detectable rates. The formation of a normal–infectious heteromeric complex lowers the activation energy barrier to convert the normally folded protein into an infectious species. This conversion leads to further recruitment of normally folded proteins in an autocatalytic process. In essence, a single infectious prion protein in a population of normally folded ones is enough to convert the whole population via autocatalytic structure replication into an all-infectious protein population, which self-assembles into long fibers.

Prions represent a form of conformational self-replication that is so far not observed in the context of RNA biology. Minimal self-switching RNA prions can serve as a model for (i) a new class of riboswitches that provide exponential feedback, or (ii) self-induced nano-units that assemble or disassemble upon stimulation. Importantly, the computational design aspect allows for context-sensitive optimization, which is necessary for applications in synthetic biology and nanotechnology. Previous research in the field of bottom-up synthetic biology has already introduced designs of small self-assembling (for a recent review see [1]), self-replicating [35], and self-polymerizing [33] systems. While all of these results are valuable in showing that such designs are indeed possible, many designs are still done by intuition and hardly adjustable to a context-sensitive implementation such as prospective artificial life forms. Along those lines we recently submitted an experimentally verified computational design of multiple self-polymerizing ribozymes in order to study the dynamics of self-interactive systems [32].

This contribution was motivated by the question of whether RNA molecules can be designed in silico to exhibit the aforementioned prionlike behavior. We show that it is indeed possible to design such an RNA prion; whether the suggested sequence really shows the exponential refolding characteristics awaits experimental verification. The RNA prion presented here is a 49-nt-long, circular RNA, designed as a bistable molecule. It thermodynamically favors one structure (S₁) if present as a monomer and the other structure (S₂) if present as a dimer.

RNA secondary structure is a good approximation for the real RNA structure because, in contrast to proteins, RNA secondary structure captures the majority of the folding free energy. Furthermore, a well-established energy model for RNA secondary structures exists that is extensively parameterized via melting experiments [39]. On the computational side, algorithms have been developed for the thermodynamic and kinetic characterization of RNA secondary structures [28]. In particular, the energetically optimal structure as well as properties of the structural ensemble at thermodynamic equilibrium can be computed efficiently. The topology of the discrete folding landscape [12], which has a strong influence on the folding kinetics, can be analyzed in detail. Several approaches model the folding kinetics of RNA secondary structure as a stochastic process with different resolution of the folding landscape [9, 45, 26] (for a review see [10]). Recently these approaches have been extended to operate on folding landscapes that change with time, as in the case of folding during transcription [19].

Let the set Ω of RNA secondary structures be restricted to those that are formed from nested isosteric base pairs (GC, AU, GU), which have a minimal hairpin loop size of 3 nt and a maximum interior loop size of 30 nt. These restrictions are broad enough to include the vast majority of known pseudoknot-free RNA conformations, and they define a set of structures for which an experimentally determined energy model E exists [29]. Most importantly from the computational perspective, these definitions allow us to compute the structure of minimum free energy (MFE) and the equilibrium partition function (Z) in O(n³) time, where n is the sequence length.

Finding the best folding path has been shown to be a NP-hard problem [30]. However, there exist fast heuristics to compute direct (shortest) paths between two structures, such as the findpath [11] method available in the ViennaRNA package, as well as indirect paths using, for example, RNAtabupath [7]. For sequences shorter than about 100 nt, paths with a minimal energy barrier can be computed exactly with RNAsubopt [46] and barriers [12]. The first program computes a list of all RNA secondary structures within a certain energy range; barriers processes a sorted list of these conformations to compute all local minima and the saddle points connecting them by a flooding algorithm.

RNA prions are bistable molecules that favor one structure (S₁) if present as a monomer and the other structure (S₂) upon dimerization. To avoid spontaneous refolding between S₂ and S₁, the energy barrier β_S2→S1 has to be very high (see Figure 1).

In our design, S₂ forms two 10-nt hairpins that are prone to hybridize via a kissing interaction as reported for the HIV-DIS loop [8], a highly conserved stem–loop sequence found in many retroviruses. Importantly, S₂ should not only stabilize other molecules in S₂ conformation at high concentrations, but actively lower the energy barrier to refold S₁ into S₂ (see Figure 1).

If these landscape properties are fulfilled, it ensures that an initial population of only S₁ molecules will not refold into S₂ unless the refolding is triggered by an external mechanism. However, as soon as a small population of molecules in S₂ conformation is present, they can catalyze the refolding of S₁ molecules into S₂.

RNA molecules with complex energy landscapes can be designed in a two-step process. First, a fast heuristic is used to generate and select promising candidates according to a cost function that specifies thermodynamic aspects of the design objective. Second, minimal refolding barriers are computed for the set of candidate molecules and contribute to the final ranking of molecules (see Figure 2 for the prion design pipeline).

RNA kissing interactions go beyond normal pseudoknot-free secondary structures as defined above, since they comprise non-nested base pairs, and their energies can therefore not be computed from the standard energy model. While we expect the energies of the intermolecular helix to be well described by standard energy parameters, it is less predictable how the (mostly entropic) contribution of the hairpin loops involved in the kiss will change. Thermodynamic stabilities of kissing interactions similar to the HIV–DIS loop, which features a 6-nt intermolecular interaction with a 9-nt hairpin loop, were studied in detail by Weixlbaumer et al. [42]. By varying the loop sequence they demonstrated that normal stacking energies can indeed be used for the intermolecular helix. More importantly, they found kissing interactions to be surprisingly more stable than other hybridization structures. Based on the measurements in [42], we compute the energy of kissing hairpins as the energy of the intermolecular helix without energy bonuses for single-base stacking of adjacent unpaired bases (dangling ends) and add a loop energy of −4.2 kcal/mol. In contrast, other intermolecular hybridization structures are penalized by the usual intermolecular initiation energy of +4.1 kcal/mol, also confirmed by Weixlbaumer et al. [42].

These energy parameters enable us to accurately evaluate the energy of S₂ • S₂, but, unfortunately, they are not sufficient for creating a consistent energy model for all intermediate conformations along the refolding paths P_{S1+S2→S1·S2} and P_{S1·S2→S2·S2}. We therefore approximate the energy barrier β_{S1+S2→S1·S2} with the findpath heuristic and compute the best energy barrier β_{S1·S2→S2·S2}, using the RNAsubopt–barriers approach with two different energy models that serve as an upper and a lower bound (see Figure 4).

The barrier for the initiation of the intermolecular hybridization β_{S1+S2→S1·S2} is computed by first finding the best path for opening the competing helix in S₁ and second computing the energy barrier for the intermolecular hybridization reaction using the standard penalty of +4.1 kcal/mol. The bonus energy of −4.2 kcal/mol is added once the full kissing region is formed, to be consistent with the energy model E_u. The energy barrier for the initiation of the duplex formation is therefore usually either the last base pair in the process of unfolding the competing helix, or the first inserted base pair towards the kissing interaction.

switch.pl [11] of the ViennaRNA package [28] was used to design bistable molecules and was modified to support the novel cost function composed of Equations 6 and 7 as well as the folding of circular RNAs. The algorithm first builds a dependence graph in order to efficiently and fairly sample RNA sequences that are compatible with both structural constraints. Since switch.pl can only design bistable and not tristable sequences, the structural constraint for the kissing interaction was specified indirectly as a sequence constraint. This reduces the number of candidate molecules, but ensures that they always have an experimentally validated, stable kissing interaction.

The chosen sequences for the kissing interaction have a similar free hybridization energy to that of the best kissing interaction examples shown in Weixlbaumer et al. [42], but differ by point mutations to be compatible with structural constraints for S₁ and S₂. Importantly, (i) the kissing hairpins cannot form intramolecular base pairs that would compete with the formation of an intermolecular kissing interaction, and (ii) the two complementary regions forming the kissing interaction should not form an intramolecular helix in S₁, which would make them inaccessible for intramolecular interactions. The asymmetric design shown in Figure 1 allows S₂ to open a shorter helical region that has a less stable free energy than the subsequent formation of the kissing interaction.

Alternatively, one could use RNAdesign [21], which can build dependence graphs for multiple structural constraints and therefore increases the search space, allowing for novel kissing interactions. This more general design attempt will be implemented upon experimental feedback for the design presented here.

The candidate molecules computed by switch.pl are subsequently ranked by the difference of barrier heights for single-molecule refolding and kissing-dimer refolding. Computation for monomer refolding is straightforward, by computing the suboptimal structures for the monomer using RNAsubopt followed by evaluation of the minimal barrier height β_S2→S1 with barriers. For the kissing-dimer interaction, suboptimal structures were computed for both energy models described above. When using the energy model E_u (see Equation 9), suboptimals were computed with the constraint that the region involved in the kissing interaction of S₁ (red in Figure 1) is unpaired (i.e., involved in the kissing interaction); according to the energy model E_l (see Equation 10), suboptimal structures were computed for the hybrid, with the constraint that the kissing region is paired, while the part corresponding to molecule 2 was kept constant in conformation S₂. The barrier for the initiation of the intramolecular hybridization is computed using findpath for opening the competing helix in S₁ and subsequently adding the duplex energies for formation of the kissing interaction.

Out of all possible sequences that are compatible with the sequence and structure constraints shown in Figure 2, 158 different sequences were returned from switch.pl using 1000 individual runs with each 2 × 10⁶ optimization steps. After postprocessing with RNAsubopt, barriers, and findpath to compute the final ranking of the generated designs, 69 sequences either turned out not to fold exactly into the ground state structures specified as input for switch.pl, or had a higher barrier for dimer refolding than for monomer refolding and were therefore excluded from further analysis. In the remaining pool of 89 sequences, 23 showed differences between monomer and dimer refolding barriers higher than 4 kcal/mol according to E_u (see Equation 9) for dimers. These sequences were visually inspected, and we selected a candidate (ACCUGGGAACCGGCGACCCAGGUUUUCGGAACCAACGUCGGAGGUUCCU) for demonstration of prion behavior, that has (i) a very high refolding barrier for monomer refolding with +16.70 kcal/mol, and (ii) an energy landscape with as few as possible competing local minima to S₁ and S₂. The dimer refolding barriers are +11.60 and +8.60 kcal/mol for E_u and E_l, respectively. In comparison, the molecule presented in Badelt et al. [3] has the same barrier for monomer refolding (+16.70 kcal/mol), but +13.60 and +9.80 kcal/mol for prion-induced refolding.

Figure 6 shows the energy profiles of the best refolding paths between S₁ and S₂, either for a single RNA monomer or for an RNA engaged in kissing interaction with another molecule. Since S₁ is the thermodynamically favored state in monomers, we show the refolding path from S₂ (−13.70 kcal/mol) to S₁ (−16.00 kcal/mol) in the top panel. The barrier of this refolding path is 16.70 kcal/mol, making a non-induced switching of conformations unlikely.

The bottom panel of Figure 6 shows the energy profile for a scenario where an intermolecular interaction is first formed between one molecule in conformation S₁ and a second in S₂, followed by intramolecular refolding of the first molecule from S₁ into S₂. S₂ is now the favored conformation, since it is stabilized by the kissing interaction. In contrast, S₁ is destabilized, since one helix cannot be formed together with the intermolecular duplex. Theoretically, there would be a second possible duplex interaction that required S₁ to open eight base pairs in two helices, but since this interaction is not thermodynamically favored, it is not depicted in Figure 6.

For the initiation of the kissing interaction, all competing intramolecular base pairs of S₁ have to open first, and then the intermolecular base pairs can form. The energy barrier for this interaction is 8.90 kcal/mol and leads to a new local minimum conformation at −34.50 kcal/mol.

For the intramolecular refolding from S₁ to S₂, we compare the two energy models discussed above. The energy model E_u returns a barrier of 11.60 kcal/mol, while the energy model E_l results in a barrier of 8.60 kcal/mol. Note also that the folding path itself is different, due to the different modeling of involved loop regions.

In this contribution we have shown that the computational design of RNA molecules that exhibit prionlike behavior is feasible, and that the computational machinery is developed enough for a rigorous analysis of the behavior of the resulting sequences. As in the original prion system, the misfolded conformation forms, via a kissing interaction, a heterodimeric complex with the native conformation. This interaction destabilizes the native conformation and triggers refolding into the misfolded conformation. Hence, we demonstrated, at least in silico, that RNA molecules possess the necessary structural capabilities for conformational replication. The calculations show that the kissing interaction drastically lowers the activation energy for refolding. Furthermore, the misfolded conformation can oligomerize. In principle, the oligomerization could inhibit the exponential growth characteristics of the misfolded conformation. Experimental results from fiber-forming dynamic combinatorial libraries [4] show that mechanical forces lead to fiber breaking, restoring the exponential growth characteristics. One difficulty in our computational design is the lack of energy parameters for complex interaction structures that occur as folding intermediates. The design process is, however, flexible and can incorporate feedback from wet lab results. We therefore envision that practical RNA designs should be refined in a few rounds of wet lab testing and adaptation of the computational models.

Conformational replication constitutes a novel regulatory mechanism possessing highly nonlinear dynamic characteristics. This type of behavior is necessary for the construction of signal-enhancing molecular circuits. Such amplifying devices are usually hard to construct with RNA. Our design could easily be coupled with interaction sites for an external signal molecule that triggers the initial refolding event. Such a device could detect a single molecular event and translate it into a large, easily detected signal. In a prospective artificial metabolism, the infectious conformation may be used to trigger transcriptional or translational riboswitches, and, if it were regulating its own transcription, the self-switching mechanism could be used for time-delayed feedback loops. The aspect of self-assembly to long fibers—which initially does not seem of particular interest, since it limits exponential growth—may provide a basic model towards the design of an RNA-based cytoskeleton. However, the current challenge is the synthesis of the RNA molecule in one particular conformation in vivo, which, for example, may be approached using techniques involving the tRNA ligase.

This work was supported in part by the FWF International Programme I670, the DK RNA program FG748004, the EU-FET grant RiboNets 323987, and the COST Action CM1304 “Emergence and Evolution of Complex Chemical Systems.”