Secondary-Structure-Informed RNA Inverse Design via Relational Graph Neural Networks

1. Introduction

RNA molecule design is essential to gene therapy, drug development, biosensing, synthetic biology, and more recently, mRNA vaccines[1,2,3]. However, in deep learning explorations, protein sequence-structure analyses have enjoyed more advances than RNAs. These advances are partly due to the significant amount of available data on protein structures [4,5,6,7,8]. Apart from the paucity of data, other major factors posing challenges in RNA inverse design are RNA’s inherent flexibility, plasticity, in some cases its ability to form distinct conformations under different environmental factors, and non-unique structure-sequence mapping.

Structure-based design of RNAs capable of folding into alternative structures, known as RNA switches or riboswitches [9] is of great interest in biotechnology and therapeutics [10]. Probably representatives of an ancient regulatory mechanism [11,12] even long before proteins [13,14,15,16] riboswitches are usually located in the 5’ untranslated region (UTR) of an mRNA. They generally consist of a ligand-binding aptamer that upon binding to the ligand or metabolite of interest can allosterically alter the structure of the rest of the riboswitch which also includes an expression platform, in turn regulating transcription or translation of downstream genes in bacteria [9,17,18,19]. Currently categorized into 55 distinct classes, riboswitches enjoy a great structural and functional diversity [20,21,22] as well as the diversity of ligands that bind to them helping cells regulate a great variety of chemicals essential to all forms of life [18,23]. It is speculated that additional classes of riboswitches exist in not only bacteria but also eukaryotes [16,18]. Riboswitches are used in many therapeutic and biotechnological applications including Riboswitch-targeting antibiotics [24,25,26,27], designer gene control [28,29,30], molecular fuse [31,32], and research tools to explore fundamental biological processes [33].

RNA inverse design has been traditionally based on RNA’s 2D or secondary structure [34]. More recent approaches such as libLEARNA [35] incorporate sequence and secondary-structure constraints to to perform partial sequence design with outstanding performance applied on riboswitches. One of the more popular 2D inverse folding tools is ViennaRNA [36], which predicts the complete sequence.

RNA computational design directly from its 3D structure can be performed by the state-of-the-art physically based molecular modeling software Rosetta [37]. Several deep learning approaches to 3D RNA inverse design have been explored in the past few years that are based on Graph Neural Network (GNN) architectures. Some of these methods are the representation-based learning approach RDesign [38] and the diffusion model RiboDiffusion [39] that uses a Transformer-based sequence module. Another GNN-based model, gRNAde [40], which incorporates concepts from autoregressive GVP-GNN [41] and protein sequence design ProteinMPNN [6] has been shown to achieve a higher native sequence recovery compared to Rosetta and RDesign on a set of benchmark sequences that include riboswitches and ribozymes [42]. The multi-state architecture of gRNAde software uses multiple structures of a single sequence in its training and design. The architectural design, however, is more geared towards aggregating multiple structural information, rather than combining them, which could lead to loss of information in certain cases. An approach that is based on combining information from multiple structural states, however, is yet to be fully explored. In this work, based on gRNAde, a secondary-structure-informed relational Graph Neural Network (GNN) is proposed for RNA inverse design with an architecture that emphasizes learning from RNAs specifically having distinct alternative structures, i.e., riboswitches. The model uses geometric-vector-perception GNN [41] layers and has a similar architecture to gRNAde [40,43], with the following alterations:

The edge embedding does not include 16 sinusoidal positional encodings. In addition, edges consist of three different types: Primary-structure edge types, secondary-structure edge types, and tertiary-structure edge types. Primary-structure edge types are defined along the sequence backbone. Secondary-structure edge types are determined by an independent software [44,45] that detects canonical and non-canonical RNA base pairing. Tertiary-structure edge types here have a specific definition and are based on clustering nucleotides based on a combination of their backbone and Euclidean distances. An example of edge types is illustrated in Figure 1, which corresponds to a section of the immature 30S ribosomal subunit of Staphylococcus aureus. Details about edge determination are explained in Materials and Methods. Each encoder layer consists of three distinct GVP-GNNs, one for each edge type. Finally, the decoder takes updated node embeddings of the last layer of the encoder, each resulting from a distinct structural scenario that the sequence can fold into, pooling them to produce the final updated node embedding. The proposed relational GNN is also $O\left(3\right)$-equivariant and message-passing components are implemented via PyTorch Geometric [46]. The Materials and Methods section describes details about input embedding and model architecture.

The paper first defines model architecture in Materials and Methods. In the Results section, the performance of the relational GNN is then presented and compared to their corresponding gRNAde counterparts under different data splits and model variations. Finally, the Discussion section contains conclusions drawn from performance analyses and future direction.

2. Materials and Methods

2.1. Data Collection

RNA structure files were downloaded from the RNASolo [48] repository on August 1, 2024. There were a total of 14,889 PDB files at resolution $\le 4$Å. The cleaned dataset contained a total of 12091 structures, corresponding to 3963 unique sequences. Sequence types were: 2612 Protein-RNA Complexes, 479 Solo RNAs, 14 DNA-RNA Hybrids, and 858 unknowns. Sequence average length was $851.74\pm 1148.69$. Only one structure was available for most sequences and the number of structures per sequence followed an exponential distribution with an average value of 3 structures per sequence. In Figure 2, parts (a) and (b) show the distributions of the sequence length and the number of structures per sequence, respectively. Similarity between structures of the same sequence was assessed by computing the average pairwise Root Mean Square Deviation (RMSD) of structures belonging to that sequence. The average value was around 1.3Å. Complete data statistics and histograms can be generated via notebooks/data_stats.ipynb. A machine-learning ready set was then produced based on the downloaded PDB files. In cases of more than one structure in the PDB file, the first one was used. Structures corresponding to the same RNA sequence were then categorized together, representing a multiple-structure single-sequence set of data (see the data availability section for the location of processed files).

2.2. Data Splits

RNA sequences were split into training, validation, and test sets. The sizes of validation and test sets were 100 each. Two different strategies were used to split the data into training and test sets. The first one was denoted as seqid. In this procedure, the aim was to avoid exposing the model to more sequences during training. This was achieved by making the training set have less diverse sequences than the test set. RNA sequences were clustered into groups based on sequence identity using command CD-HIT [49] (threshold 90%) and those with less intra-sequence RMSD among their structures were favored for the training set. In the second strategy, the aim was to avoid exposing the model to more structures during training. This was achieved by making the training set have less diverse structures than the test set. RNA sequences were clustered into groups based on structural similarity using command qTMclust [50] (threshold 0.45%) and those with less intra-sequence RMSD among their structures were favored for the training set. Both data splitting procedures were according to [40].

2.3. Metrics

Performance parameters were perplexity, accuracy, recovery, and sc score. Perplexity reflects the ambiguity within nucleotide probabilities. For instance, a model that predicts equi-probable nucleotides for each node is undesirable. Lower output perplexity is generally favorable1. Accuracy is the ratio of correctly predicted nucleotides. Recovery is the ratio of correctly predicted nucleotides from a population of sequences that are sampled from output probabilities. Here, 16 predictions were used. Similarly, sc score is the secondary-structure self-consistency score averaged over 16 predictions. Software Eternafold [51] was used to measure how well the predicted secondary structures of sampled output sequences match that of the ground-truth sequence.

2.4. Graph Construction of RNA Structures

The 3D coordinates of atoms of an RNA structure, taken from its corresponding PDB file, were used to construct a geometric graph representation (nodes and edges) of the molecule. Three types of edges were defined and denoted as primary-structure ${{\overrightarrow{\mathit{e}}}^p}_{i,j}\in E^p$, secondary-structure ${{\overrightarrow{\mathit{e}}}^s}_{i,j}\in E^s$, and tertiary-structure ${{\overrightarrow{\mathit{e}}}^t}_{i,j}\in E^t$ edge types between nodes i and j. Edge types are mutually exclusive: $E^p\cap E^s=\varnothing$, $E^p\cap E^t=\varnothing$, $E^s\cap E^t=\varnothing$.

Nodes

The 3D coordinates for P, C4’, N1 (pyrimidine), or N9 (purine) atoms of each nucleotide were used to define the corresponding node for that nucleotide [52]. The C4’ coordinates of nucleotide i, ${\overrightarrow{\mathit{x}}}_i\in {\mathbb{R}}^3$, represented the 3D coordinates of the corresponding node i. The 3D coordinates of P, C4’, and N1/9 as well as the atoms of the adjacent nucleotides on the RNA backbone were used for feature representation.

Primary-Structure Edge Types ${\overrightarrow{\mathit{e}}}_{i,j}^p$

Pairs of consecutive nucleotides on the sequence, i.e., $j=i+1$, were assigned primary-structure edge types in both forward ${\overrightarrow{\mathit{x}}}_{i+1}-{\overrightarrow{\mathit{x}}}_i$ and backward ${\overrightarrow{\mathit{x}}}_i-{\overrightarrow{\mathit{x}}}_{i-1}$ directions.

Secondary-Structure Edge Types ${\overrightarrow{\mathit{e}}}_{i,j}^s$

Pairs of nucleotides in secondary-structure RNA base pairings, $(i,j)\in BP$, were assigned secondary-structure edge types. For each base pairing $(i,j)$, two edges ${\overrightarrow{\mathit{x}}}_i-{\overrightarrow{\mathit{x}}}_j$ and ${\overrightarrow{\mathit{x}}}_j-{\overrightarrow{\mathit{x}}}_i$ were defined. Here, secondary structure base pairings included both canonical and non-canonical pairings [53] and determined by feeding the corresponding PDB file of an RNA to X3dna-dssr [45].

Tertiary-Structure Edge Types ${\overrightarrow{\mathit{e}}}_{i,j}^t$

Pairs of nucleotides that did not form secondary structure base pairings but were in close proximity of each other were considered for having what is referred to here as tertiary-structure edge types. The procedure was as follows: Given the secondary-structure base pairs, the unpaired nucleotides were first identified, i.e., $i,j\notin BP$, and then clustered using Density-Based Spatial Clustering of Applications with Noise (DBSCAN). Default parameters for DBSCAN were DBSCAN_eps = 20Å and min_samples = 5. Subsequently, those pairs of nucleotides that belonged to the same cluster with distance of $d={∥{\overrightarrow{\mathit{x}}}_i-{\overrightarrow{\mathit{x}}}_j∥}_2<2\times$DBSCAN_eps, and were also far apart on the RNA sequence were selected. Having distance $|j-i|\ge 500$, or at least 500nt, was considered far apart. i.e., default primary_dist = 500. The number of original tertiary-structure edge types increased exponentially with sequence length. The edges were randomly selected from the original sample, such that the number of tertiary-structure edges was strictly less than or equal to $4\times L$, where L is the length of the corresponding RNA sequence.

The number of edges selected from each RNA molecule is shown in Figure 3). As we can see, the number of edges selected is linear with respect to the length of the sequence. The dashed magenta line shows the cut-off value for cases where the number of original edges was more than the threshold. Edges of RNAs with higher counts were down-sampled to strictly match the threshold.

2.5. Relational Graph Neural Network

Input Embedding

Based originally on [41,54] and identical to the implementation presented in [43], node features were derived from three distances, three inter-atom angles, and three dihedral angles. The number of scalar features associated with each node i was then: sin and cos values of the six angles, concatenated by the three lengths, totaling 15 values. The vector features associated with each node i were four unit vectors corresponding to forward and backward C4’ - P, and forward and backward C4’ - N1/9 vectors. Scalar and vector features of node i are denoted as $\left({\mathit{s}}_i,{\overrightarrow{\mathit{v}}}_i\right)$. Edges of all types, ${\overrightarrow{\mathit{e}}}_{i,j}$, had identical features. Edge features consisted of both unit vectors and scalar values. Normalized vectors $\langle {\overrightarrow{\mathit{e}}}_{i,j}\rangle ={\overrightarrow{\mathit{e}}}_{i,j}/∥{\overrightarrow{\mathit{e}}}_{i,j}∥$ were used as unit vectors. Scalar features consisted of ${∥{\overrightarrow{\mathit{x}}}_i-{\overrightarrow{\mathit{x}}}_j∥}_2$ encoded by 16 RBFs. Let ${\mathit{x}}_i\in X$ represent both scalar and vector features of node i and ${\mathit{e}}_{i,j}\in E$ represent both scalar and vector features of edge $(i,j)$. Hence, $E=\{E^p,E^s,E^t\}$ Let matrix T represent the type of each edge in E in its corresponding indices. For a single sequence-structure scenario, each RNA structure is then represented as a graph $G=(X,E,T)$. In case of multiple (K) conformations, the graph corresponding to structure k is denoted as $G^{\left(k\right)}=(X^{\left(k\right)},E^{\left(k\right)},T^{\left(k\right)})$. All structures belonging to the same sequence, $\{G^{\left(1\right)},\cdots ,G^{\left(K\right)}\}$, are then passed into the GNN model in a single batch. Unlike [40], however, where the K graphs were merged into a single multi-graph with total edges being union of those of individual structures, here edges of each structure k are applied to the model separately. Hence, although the number of nodes is identical amongst the K graphs, the number of edges is not. The rationale behind the above choice of input representation of multiple structures was to favor the expressiveness of an RNA molecule with drastically different conformational states, e.g., a riboswitch, over that of an RNA molecule with only one innate structure but having multiple experimental data on that structure.

Model Encoder

The overall architecture of the pipeline consisted of a linear model, L stacked layers of relational GNNs, and followed by a decoder layer which produces nucleotide probabilities for each sequence location i. Node and edge embeddings are first passed through the linear model. Encoded features are then passed to the first relational GNN layer. The relational GNN consists of three $O\left(3\right)$-equivariant GVP-GNN [41] components. Using edge-type indices, T, the encoded features of $G^{\left(k\right)}$ are first broken down to three sub-graphs, ${G^p}^{\left(k\right)}=(X^{\left(k\right)},{E^p}^{\left(k\right)})$, ${G^s}^{\left(k\right)}=(X^{\left(k\right)},{E^s}^{\left(k\right)})$, and ${G^t}^{\left(k\right)}=(X^{\left(k\right)},{E^t}^{\left(k\right)})$, each representing a graph of a specific edge type. Each of the sub-graph encodings is fed into a separate GVP-GNN (eq. 1), and then pooled using learned weights (eq. 2) to produce updated node features (eq. 3):

$$\begin{array}{cc}\hfill {{\mathit{m}}_i}^p,{\overrightarrow{\mathit{m}}}_i^p& :=\sum _{j\in {\mathcal{N}}_i^p}MSG\left(\left({\mathit{s}}_i,{\overrightarrow{\mathit{v}}}_i\right),\left({\mathit{s}}_j,{\overrightarrow{\mathit{v}}}_j\right),{\mathit{e}}_{ij}^p\right)\hfill \\ \hfill {{\mathit{m}}_i}^s,{\overrightarrow{\mathit{m}}}_i^s& :=\sum _{j\in {\mathcal{N}}_i^s}MSG\left(\left({\mathit{s}}_i,{\overrightarrow{\mathit{v}}}_i\right),\left({\mathit{s}}_j,{\overrightarrow{\mathit{v}}}_j\right),{\mathit{e}}_{ij}^s\right)\hfill \\ \hfill {{\mathit{m}}_i}^t,{\overrightarrow{\mathit{m}}}_i^t& :=\sum _{j\in {\mathcal{N}}_i^t}MSG\left(\left({\mathit{s}}_i,{\overrightarrow{\mathit{v}}}_i\right),\left({\mathit{s}}_j,{\overrightarrow{\mathit{v}}}_j\right),{\mathit{e}}_{ij}^t\right)\hfill \end{array}$$

(1)

$${\mathit{m}}_i,{\overrightarrow{\mathit{m}}}_i:=\sum _{e\in \{p,s,t\}}{w}^e\left({{\mathit{m}}_i}^e,{\overrightarrow{\mathit{m}}}_i^e\right)$$

(2)

$${\mathit{s}}_i^{\prime },{\overrightarrow{\mathit{v}}}_i^{\prime }:=UPD\left(\left({\mathit{s}}_i,{\overrightarrow{\mathit{v}}}_i\right),\left({\mathit{m}}_i,{\overrightarrow{\mathit{m}}}_i\right)\right)$$

(3)

Messages with p, s, and t indices reflect primary-structure, secondary-structure, and tertiary-structure messages from neighboring nodes, respectively. Learned weights, $\{w^p,w^s,w^t\}$ reflect the significance of each edge type, i.e., the attention given to different edge types. The above message-passing and pooling mechanism applies to any relational GNN layer and layer indices are omitted for simplicity. For instance, edge-type weights of layer l can be denoted as: $\{{w^p}^{\left(l\right)},{w^s}^{\left(l\right)},{w^t}^{\left(l\right)}\}$. In order to assess the impact of the secondary and tertiary edge types separately, two different versions of the model were used, one including primary and secondary edge types only, denoted by suffix -2D, and one including all edge types, denoted here by suffix -3D.

Scalar and vector node features are updated at every layer (eq. 3). For instance, the output of layer l is $\left({\mathit{s}}_i^{\left(l\right)},{\overrightarrow{\mathit{v}}}_i^{\left(l\right)}\right)$. The above message passing scheme is applied to all K structural conformations in the input. At the end of the last encoder layer, the updated node embedding for the k-th structure is $\left({\mathit{s}}_{i,k}^{\left(L\right)},{\overrightarrow{\mathit{v}}}_{i,k}^{\left(L\right)}\right)$. The final node (or nucleotide) embedding given all K structures of the RNA sequence is then an average of individual estimates:

$${\mathit{s}}_i,{\overrightarrow{\mathit{v}}}_i:=\sum _{k\in K}\left({\mathit{s}}_{i,k}^{\left(L\right)},{\overrightarrow{\mathit{v}}}_{i,k}^{\left(L\right)}\right)$$

(4)

Model Decoder

The original model [40] presented two different decoder implementations: Non-AutoRegressive, NAR-original, and autoregressive AR-original. The NAR-original decoder consisted of two GVP-GNN modules and decoded the node embedding into a dimension of 4 to reflect the four nucleotide probability estimates at each sequence position. The AR-original decoder was an autoregressive model [55] that updated each node embedding based on the nodes there preceded it in the sequence. Here, based on the above two decoders, two secondary-structure-informed decoders by the names NAR-informed and AR-informed were deployed. The implementation of NAR-informed was identical to NAR-original. In the AR-informed decoder, node embeddings specific to structure k are first independently updated, $\left({\mathit{s}}_{i,k},{\overrightarrow{\mathit{v}}}_{i,k}\right)$. A final decode layer then pools the individual embeddings to generate final nucleotide probabilities for each sequence position.

3. Results

The proposed RNA inverse design is mainly based on defining different edge types. Unlike primary- and secondary-structure edge types, however, the tertiary-structure edge types are not as deterministic and vary based on the choice of modeling. Here, the intention behind tertiary-structure edges was to capture some of RNA molecule’s tertiary patterns that are evolved beyond the secondary structural constraints. Here, a simple model is used based on the relative orientation of unpaired nucleotides. A general analysis of all the tertiary-structure edge types is provided in  A.1. Given the unevenness of the relative-orientation distribution of the unpaired nucleotides (Figure A1 part (a)), it was speculated that introducing tertiary-structure edge types as defined above may contain certain hidden information about RNA nucleotide positioning in 3D space.

The secondary-structure-informed relational GNN architecture use is shown in Figure 4. The model consists of L relational message-passing GNN layers. Each layer pools messages from three distinct GVP-GNN components each corresponding to an edge type. Note that a single GVP-GNN component can be used instead of three where, ${\mathit{e}}_{ij}$ would be defined as concatenation of ${\mathit{e}}_{ij}^p$, ${\mathit{e}}_{ij}^s$, and ${\mathit{e}}_{ij}^t$. In fact, Eqs (1) and (2) are mathematically equivalent to Equation (1) of [40]. However, since edges are mutually exclusive, the merging of edge types would result in a sparse adjacency matrix. In addition, each ${\mathit{e}}_{ij}$ vector would have been 3 times an original edge vector, padded by zeros in corresponding locations. Although a single GVP-GNN would have been equivalent to the three original GVP-GNNs, deploying three GVP-GNNs has the advantage of requiring a smaller input size for edges.

Given multiple structures of a particular sequence, each individual sequence-structure pair is fed into the encoder, separately, resulting in a k set of final updated node embeddings (output of layer L). The node embeddings, each resulting from a distinct structure, are then pooled before being fed into the decoder, which produces the final probabilities $p\left(i\right)$ for each nucleotide i:

$$p\left(i\right)=(1/K)\sum _{k\in K}p\left(i\right|structur{e}_k)$$

(5)

Apart from the above re-structuring, The GVP-GNN components, linear, and final decoder layers are adopted from and according to the gRNAde software [40,43]. Two different decoders were introduced in [40] by names one-shot and autoregressive, which are referred to here as NAR-original and AR-original, respectively. Implementations for both models were altered to implement the above secondary-structure-informed relational GNN versions, denoted as NAR-informed and AR-informed, respectively. For each altered implementation, two different versions of *-informed-2D and *-informed-3D were considered. In the case of AR-informed-*, the pooling of updated node embeddings was performed after autoregressive node prediction of individual structure. Hence, eq. (5) holds for all implementations.

Table 1 contains a performance comparison between the original gRNAde and the secondary-structure-informed relational graph neural network under different data split strategies and for different parameters such as the number of encoder layers and inclusion/exclusion of tertiary-structure edge types.

In addition to the [45] software, other choices of software such as Fr3d [56], and corresponding implementation NA_pairwise_interactions.py were considered, which led to similar results (data not shown). The Fr3d software annotates base pair interactions in RNA 3D structures according to the Leontis–Westhof classification [53].

Performance measures of Table 1 were obtained from single runs. Therefore, variation in performance is not available. Based on table (1) of [40], however, we can see that the upper bound of the standard deviation of mean Recovery for $k=3$ is $\pm 1.19$. Since models here are very similar to [40], we assume the mean Recovery to have similar standard deviation measures.

The original autoregressive model AR-original trained on seqid data split, and implemented as part of gRNAde had a lower perplexity as well as higher sc score across all models. Higher sc score is consistent across both sequence-based seqid and tertiary-structure-based structsim data splits.

Comparing models based on accuracy and recovery metrics, however, shows a higher performance for the *-informed-* models compared to their *-original counterparts. The higher performance is true regardless of the choice of data split strategies. In fact, the *-informed-* models only have 2 layers while the original model is set to 4 layers. The comparison here is based on the assumption that the mean Recovery deviation is bounded by $\pm 1.19$. Finally, increasing layers from 2 to 4 further improves the performance of the AR-informed-3D model to recovery=58.37.

Comparing *-informed-2D model performances to their *-informed-3D counterparts shows that tertiary-structure edge types did not result in any major improvement.

Inverse design performance of the seqid-trained AR-informed-3D model was also evaluated on the 14 RNA structures of interest identified by ([42], Supplementary Table 2) and compared to previous methods (See Table 2). Recovery values are based on averaging prediction accuracy over 16 predicted sequences in each case. Our proposed secondary-structure-informed approach had an average of 0.6044 native sequence recovery rate, slightly higher than gRNAde (0.5682) and also higher than other methods including Rosetta (0.45) and RDesign (0.4296). All Recovery values are according to [43].

In order to showcase an example of designing an RNA, an experimentally derived structure of thiamine pyrophosphate (TPP)-sensing riboswitch was selected and passed into the structsim-trained AR-informed-3D model to obtain 16 sequence predictions. The crystal structure of the TPP riboswitch, PDB ID 8F4O, was resolved in the ligand-free (apo) state and with resolution 3.1Å [57]. The sequence with the highest recovery rate (here, 0.7692) was then selected as the final prediction, referred to as Pred_1. Table 3 shows the true and predicted sequences. The structure of Pred_1 was then computationally determined via AlphaFold 3 [58]. AlphaFold provided 5 different predicted structures, the first of which had the highest similarity (lowest RMSD) with the true structure of the TPP riboswitch aptamer. Figure 5 shows the true and predicted structures, in orange and green, respectively. The first and last nucleotides of both structures are shown in licorice mode for ease of visualization. Aligned structures had RMSD of (15.88Å) among all their atoms and RMSD of (16.38Å) among their C4’ atoms. The disconnects in the true structure are due to numbering schemes in the corresponding PDB file. Although the two structures do not have perfect alignment, Figure 5 shows an extended Y-shaped conformation for the predicted structure (green), similar to a typical TPP riboswitch.

4. Discussion and Conclusions

In this work, a geometric GNN-based model was presented to improve the computational design of RNA, especially those with multiple structural states, such as riboswitches. The proposed model named secondary-structure-informed RNA inverse design is mostly taken from gRNAde, with some major differences. These differences were the use of different edge types to represent RNA structural graphs and explicitly combining results from multiple structures in corresponding sequence prediction. Details about these differences are explained in Materials and Methods. It was shown that the proposed model does indeed improve both accuracy and native sequence recovery compared to gRNAde. Using the seqid training/test sets, the autoregressive AR-informed-2D with only two layers obtained 71.12% sequence recovery, while its gRNAde counterpart with four layers produced 66.08% recovery. Similar improvements were seen under the structsim training/test sets (comparing 68.61% with 62.72%) as well as over different choices of decoders (see Table 1). Results demonstrate that the addition of secondary-structure base-pairing information in RNA inverse design and combining predictions via a relational graph architecture does indeed improve performance since both models share similar architectural components and parameters. Improvements may also be due to the fact that identification of canonical and non-canonical base pairings, namely the 12 possible orientations [53], are done independently, relieving the GNN model from the burden of identifying such edges.

Performance of AR-informed-3D was also evaluated on a different benchmark and compared to different methods. The dataset presented in [42] contains 14 structures which include riboswitches and ribozymes. It is shown in Table 2 that the proposed model achieves higher average native sequence recovery on this set over previous methods. The AR-informed model achieved 60.44 native sequence recovery rate, higher than gRNAde (56.82), as well as Rosetta (45) and RDesign (42.96). These results highlight the ability of the proposed approach in the inverse design of novel structures.

4.1. Future Work

In order to identify RNA base pairings, we used X3dna-dssr [45]. Therefore, unlike gRNAde which is a stand-alone software, our model relies on an external component to identify secondary structural base pairing patterns. Performance measures presented here are only valid if the secondary structure of the desired RNA can be independently determined. The gRNAde has the advantage of being able to handle the general case of only using the 3D structure to design sequences. In order to eliminate the need for external software, Python packages such as rnaglib [59] can be used which provide a graph representation of RNA 3D structures with edges representing base pairs and spatial interactions.

Metrics recovery and sc score evaluate primary- and secondary-structure similarity with the ground truth and are only generic assessments. Two RNAs with dissimilar sequence and secondary structure predictions might still have a very similar tertiary structure. A gold standard metric for an effective RNA inverse design should be based on similarity in key tertiary structural information, which was partly addressed by [40]. Reliable 3D structural prediction tools such as the state-of-the-art tools [58] combined with effective similarity measures such as RMSD or novel metrics inspired by the protein realm [60], can aid us in such comparisons. For the case of RNAs with more than one structural state, however, a more complex metric may be needed that reliably assesses the sequence’s ability to fold into the two desired alternative structures of a typical riboswitch.

The autoregressive AR-informed model led to the highest performance, supporting the fact that including auto-regression is very effective in RNA inverse design. In this work, auto-regression and pooling are done for each structure separately, only to combine predictions at the last step eq. (5). One direction for improvement would be to have more effective auto-regressive techniques. for instance, we can have auto-regression and pooling combined for each iteration i, i.e., determining $p\left(i\right)$ from all K structures and only then repeating the process for $p(i+1)$. Exploring novel, possibly non-causal, auto-regressive strategies may also improve the inverse design of RNAs with multiple structures.

The incorporation of tertiary-structure edge types here did not lead to any major performance increase. This may be due to the fact that our tertiary edge determination was too simplistic and could not capture any meaningful tertiary dependencies in RNA structures. Indeed, RNA 3D structural data such as base pair geometry and protein binding sites have been successfully used to predict the corresponding RNA function [59], which highlights the usefulness of extracting more complex structural patterns from RNAs. More effective modeling of RNA complexes, for example, ligand-binding configurations, may provide a clearer answer to whether learnable and generalizable patterns can be extracted by machine learning that could improve more intricate RNA inverse design challenges.

Appendix A.1. An Exploration into Tertiary-Structure Edge Type Characteristics

Using the all RNA structures available in the dataset, all nucleotide pairs having tertiary-structure edge types between them were identified. Euclidean distance and relative orientation between each of the corresponding node pairs were considered for further analysis. Value d was a scalar representing Euclidean distance between the C4’ coordinates of the two nucleotides. Filtering for Euclidean distance of $2\le d\le 20$ in Åresulted in around 43 thousand such pairs. As illustrated in Figure A1(a), the relative orientation was defined in terms of spherical coordinates with $\psi$ being the angle from the positive z-axis and $\theta$ being the angle between the radial line and a polar axis. Angle $\psi$ was defined as the angle between normal vectors of the two (C4’,P,N1/9) planes of the two nucleotides. Angle $\theta$ was defined as the angle between the two C4’-P vectors. The density plot of the normalized relative orientation of the nucleotide pairs is shown in Figure A1(b). Red represents high-density and blue represents low-density orientations. Certain relative pairwise orientations were extremely rare in 3D space. In particular, relative orientations with angle $\theta \approx 90$ were not found, regardless of values of $\psi$ or d. On the other hand, certain orientations such as $\psi =70$ were more prevalent representing an uneven distribution of relative orientations.

Figure A1. Relative orientations of distant nucleotide pairs having tertiary-structure edge types. a) Angle $\psi$ was defined as the angle between normal vectors of the two (C4’,P,N1/9) planes of the two nucleotides. Angle $\theta$ was defined as the angle between the two C4’-P vectors. b) Density plot of the normalized relative orientation of around 43 thousand nucleotide pairs, taken from RNA structures of the old Dataset. All pairs had Euclidean distance of $2\le d\le 20$ in Å. Red represents high-density and blue represents low-density orientations.

Figure A1. Relative orientations of distant nucleotide pairs having tertiary-structure edge types. a) Angle $\psi$ was defined as the angle between normal vectors of the two (C4’,P,N1/9) planes of the two nucleotides. Angle $\theta$ was defined as the angle between the two C4’-P vectors. b) Density plot of the normalized relative orientation of around 43 thousand nucleotide pairs, taken from RNA structures of the old Dataset. All pairs had Euclidean distance of $2\le d\le 20$ in Å. Red represents high-density and blue represents low-density orientations.