Qualitative Modeling of RNA Structure

1 Introduction

RNA is a single-stranded molecule consisting of a se­ quence of bases. However, the underlying linear struc­ ture of RNA is embedded in a complex secondary struc­ ture because in three-dimensional space it is possible for bases in the molecule to form chemical bonds with other bases in the molecule which are not neighbors in the underlying linear structure. Determining the secondary structure -("folding") of an RNA molecule from its un­ derlying linear structure is a complex problem involving both spatial reasoning and the use of knowledge of chem­ istry and biology.

Both spatial reasoning and scientific modeling have been studied before by artificial intelligence researchers. However, the problem of determining RNA secondary structure does not fit smoothly into the traditional qual­ itative physics framework. We believe that this domain presents challenges which require new research in both spatial reasoning and reasoning about physical systems.

RNA forms a folding structure by base pairing be­ tween complernentary bases of the molecule. The struc­ ture can be determined directly by X-ray crystallogra- An RNA molecule is made up of ribonucleotides linked together in a chain. Each ribonucleotide comprises one of four different bases, adenine (A), cytosine (c), gua-­ nine (G), and uridine (u) . Single-stranded RNA often folds back onto itself in structures stabilized by hydro­ gen bonds between A and u and between G and c. The A-U and G-C base pairs are called canonical pairs. The so­ called G-U wobble pair is also frequently found although it may have a weaker bond than the canonical pairs. These two types of parings are considered between any two bases and all the other types are considered as mis­ matches.

We make the following assumptions to solve the prob­ lem:

1. At least two homologous sequences are available and they are properly aligned.
2. Homologous RNA sequences have similar folding structures.

The folding structure is described in a discrete model because continuously varying parameters such as bond lengths or angles are not considered; either a hydrogen bond exists between two bases or it does not. Once the connectivity relation is determined, the shape of a secondary structure is unambiguously determined as we show in the next section.

3 Modeling the Structure Formation Process

In the RNA folding structure, each structure element contributes to the total struCtural stability. Helices formed by consecutive base pairings contribute to in­ creased stability, while single stranded parts destabilize neighboring helices. Many experiments have been car­ ried out with chemically synthesized RNAs to obtain the structural stability measured by the loss in free energy in kcaljmol [Tinoco et al., 1973]. A data set offree energy values describing the contribution of each structure ele­ ment is called an energy model. Several energy models are available but none of them can be considered per­ fect. Predicted structures are often radically different if a different energy model is used or a slight change is made to the energy values. Since the energy model is be­ lieved to be inaccurate as well as incomplete, we do not make commitment to the numeric values of the energy model. Instead, we use approximation of the model and seek other possible sources of knowledge, such as homol­ ogous RN As, biochemical experimental data, and biolog­ ical property of RNA. The knowledge from biochemical experiments or biological properties is not assumed to be available all the time, but it can be incorporated into the model whenever it is available. Having such prior knowledge enables pruning the search space.

The reasoning process of our model proceeds in three phases: (1) comparative analysis of homologous se­ quences, (2) construction of structure model, and (3) simulation of the process of structure formation over time. Reasoning process of the model is made explicit so that the sensitivity of the model to parameter values can be examined. In this section, we describe each of the reasoning process in detail.

3.1 Comparative Analysis of Homologous Sequences

Homologous sequences are analyzed to identify potential helices, which usually maintain complementarity in the homologous sequences. Suppose that an RNA sequence contains a subsequence GAGCU. GAGCU has the potential to base pair with any AGCUC within the RNA and such a subsequence may occur several times in that RNA. In order to identify which occurrence of AGcuc actually pairs with GAGCU, the homologous RNA sequences from different organisms are examined.

sequence 1: GAGCU ... AGCUC ... AGCUC ... AGCUC sequence 2: GAUCU ... AGCUC ... AGCUC ... AGAUC sequence 3: GAUCU .. . AGCUC ... AGCAC ... AGAUC sequence 4: GACCU ... AGCUC ... AGCAC ... AGGUC

In the above example, only the last occurrence of AGCUC changes to maintain the ability to pair with GAGCU. Such. "compensating base changes" at two po­ sitions are considered as evidence of the presence of a helix. Given an alignment of M homologous sequences of maximum length N , the program examines the se­ quences to identify potential base pairings based on the following algorithm:

for each base i in the alignment for each base j ( > i) in the alignment determine the relation BP( i, j).

The relation BP( i, j) represents the type of base pair­ ing between base i and base j for all the M sequences in the alignment, not for a single sequence. Since differ­ ent sequences can have different types of base parings, the type of canonical or wobble pair cannot be used for this relation. Instead, we define the relation BP( i, j) as follows:

1. If base i pairs with base j in all the sequences and there is no variation in both i and j , B P ( i, j) is an exact-invariant pair (denoted by o ).
2. If base i pairs with base j in all the sequences and there is a compensating base change in i and j, BP( i, j) is an exact-variant pair ( * ).
3. If base i forms a G-U wobble pair with base j in most sequences, BP(i,j) is a wobble pair (w).
4. If base i pairs with base j in most sequences but not in all the sequences and the frequency of the mismatches does not exceed a specified number, BP( i, j) is an inexact pair ( + ).
5. If the frequency of mismatches between base i and base j exceeds a specified number, BP( i, j) is a mis­ match (·).

17-----------------

Figure 1: (A) Alignment of 5 homologous sequences and their consensus sequence. (B) Covariation matrix of the sequences. (C) Folding structure inferred from the ma­ trix. (D) Structure prepared with a drawing program.

The result of the comparative analysis is represented in an N x N covariation matrix, where entry ( i, j) repre­ sents a relation BP(i,j). Since the relation BP(i,j) is symmetric, only upper triangular part of the matrix is used. In the covariation matrix a potential helix of RNA appears as a diagonal of non-mismatch symbols ( o, *· w, or +) from 1~pper right to lower left. Fig. 1A shows the alignment of five RNA sequences and Fig. 1B shows the resi.1lt of the comparative analysis of them. The last se­ quence in the alignment is a consensus sequence of the five homologous sequences, derived by the program. Up­ per case letter in the consensus sequence represents a base conserved in all the sequences and lower case repre­ sents a base which is not conserved in all the sequences but most representative one in the position. The con­ sensus sequence is used for drawing the structure in this example. Although there are several non-mismatch sym­ bols in the matrix of Fig. 1B, there is only one potential helix (1-5/13-17) which consists of at least 3 base pair­ ings, and the folding structure inferred from the matrix is shown in Fig. 1C, represented in a form of matching parentheses. The representation is directly readable by a drawing program for a better quality representation 2 (Fig. lD). Notice that any L x L submatrix (L :S N) of the covariation matrix contains potential helices at a point when the RNA sequence grows to have L bases. As we show below, we use a series of submatrices to model the structures at intermediate stages during the RNA sequence growth.

3.2 Construction of Structure Model and Simulation

In the previous section, we discussed the use of a co­ variation matrix in finding a folding structure. In many cases, however, a folding structure cannot be determined directly from the matrix due to conflicts among potential helices. For example, consider a helix (1- 5/14- 18) and a helix (1-5/19- 23) in the covariation matrix of Fig. 2B. The two helices conflict each other because bases 1-5 cannot pair with both bases 14-18 and bases 19- 23 at the same time (first folding constraint in Section 2) . The system uses a heuristic algorithm with spatial constraint propagator incorporated in it. The algorithm first identi­ fies all the potential helices which are not self-conflicting (i.e. the helix does not contain any base pairing ( i, i) and the length of hairpin loop is at least 3). For each potential helices identified, two types of scores are com­ puted, which are weighted sums of five parameters: helix length, number of exact-variant pairs, number of wobble pairs, number of inexact pairs, and the length of a pos­ sible hairpin loop created by the helix. sl is computed for finding the most ~table hairpin loop-stem structure of. the overall structure; and s2 is for measuring the "good­ ness" of a helix in general positions (e.g. helix next to a bulge loop, internal loop, multiple loop, or hairpin loop). A helix with the highest S1 value is called "besthelix" in the following algorithm. If there is any other helix which forms the same stem structure with the besthelix (i.e. these helices are intervened only by a bulge or an inter­ nal loop), it becomes a supporting helix of the besthelix. The best helix and its supporting helix (if any) consti­ tute a skeleton of the overall structure. Once a skeleton is found, it is used to eliminate the helices which violate the spatial constraints with respect to the skeleton. Each of the remaining helices is examined in decreasing order of s2, and is attached to currently maintained structures if it fits them. If it cannot fit, an alternative structure is created for it. The alternative structure initially con­ tains the skeleton and the helix which had conflicts. The heuristic algorithm is briefly stated below. Propagation of the folding constraints occur in step 4 and 6 of the algorithm. The knowledge about the scoring functions comes from energy model and comparative analysis of homologous RN As, which will be described in detail in section 4.

2 The RNA structures in this paper were drawn with LoopViewer by D. G. Gilbert at Indiana University.

1. bestscore = -oo
2. for each potential helix i identified compute sl and s2 if (S1 > bestscore) besthelix = helix i bestscore = sl of helix i
3. skeleton = {besthelix} U {any supporting helix}
4. Remove the helices which conflict with the skeleton
5. Sort the remaining helices with respect to S2
6. for each of the sorted helices if it fits current structures, add it to them otherwise, create an alternative structure for it

Our system finds 1 - 5 common foldings of M sequences of length N in O(MN 2 + N 3 ) time and O(N 2 ) space. This is the same requirement as an ordinary dynamic program for folding a single sequence [Zuker and Stiegler, 1981]. For the process modeling, structures at each in­ termediate stage of L bases ( L :S N) can be found in O(N 3 ) time and O(N 2 ) space. If the series of structures are obtained by examining the covariation matrix at k intermediate stages, the total time and space is O(k ·N 3 ) and O(k · N 2 ), respectively. If k is limited to a constant, process modeling does not raise the overall complexity. Detailed analysis of the complexity appears in [Han and Kim, 1993].

4 Experimental Results

All the algorithms described above have been imple­ mented in C on a Sun SparcStation running SunOS Re­ lease 4.1.1. The program is controlled by 8 parameters. Fig. 3 describes parameter values used for the experi­ ment in Fig. 2.

number of homologous sequences = 10 minimum length of a helix = 3 number of mismatches allowed = 1 weight values of a helix for S2 : length = 1.0 frequency of exact-variant pairs = 0.5 frequency of wobble pairs = -0.5 frequency of inexact pairs = -2.0 length of a possible hairpin loop = -0.01

In the covariation matrix (Fig. 2B) there are 63 di­ agonals which have at least 3 base parings. All the di­ agonals are potential helices, but they cannot coexist at the same time due to the folding constraints in Sec­ tion 2. Fig. 2C shows three possible structures of TAR inferred from the covariation matrix, which are most sta­ ble based on the scoring function. The first structure in Fig. 2C is identical to the one predicted by a biochem­ ical experiment [Cullen, 1990), which has been corrob­ orated by many other experiments later. Actually the superiority of the first structure to the other two can be determined by looking at them by "eye". In the sec­ ond structure, the difference of unpaired bases on either side of an internal loop is too large (11 - 1 = 10), and this will greatly destabilize the adjacent helices. In the third structure, the two helices ( GgU 1 ACC and ucucuG­ GUUAG/CUAaCUAGGGA) are identical to those of the first structure, but the other two helices in the third structure are less stable than those of the first structure in terms of the length of a helix and the number of wobble pairs.

Fig. 4 shows a series of TAR structures at nine in­ termediate stages. At stage 7 and stage 9, there are 2 and 3 structures possible, respectively. In general there are four possible transitions from a structure to another: (1) forming a new helix is initiated (e.g. from stage 3 to 4 and from stage 5 to 6 in Fig. 4); (2) a helix in a structure is broken down and a new helix is formed (e.g. from stage 1 to 2, from stage 6 to 7, from stage 7 to 8, and from stage 8 to 9); (3) an existing helix is elongated (e.g. from stage 2 to 3); and ( 4) additional bases are at­ tached to single stranded part (e.g. from stage 4 to 5). Transitions from stage 7 to 8 and from 8 to 9 seem to be most difficult in terms of the amount of base pairings to be broken down. Current system finds the most stable structures in each stage instead of the ones which have the easiest transition from the previous stage.

To obtain a stable common structure, it is generally better to look at more sequences. However, if there is little variation between the sequences, a large number of sequences are not necessary. The value of minimum length of a helix greatly affects the number of potential helices. The smaller the value, the more the potential helices. It also influences the prediction. If the num­ ber of minimum length of a helix is set to a larger value than 3, then a rough structure is obtained where helices shorter than the threshold value are missing but all the other helices stay. Therefore, increasing the minimum length of a helix does not affect much the prediction other than that shorter helices are ignored. Decreasing the number of minimum length of a helix from 3 to 2 has a slightly more complicated influence on the prediction; it can simply add new helices of length 2 to the previous structure while preserving all the existing helices in it; helices of length 2 can sometimes, however, defeat longer helices in the previous structure and replace them. The value of number of mismatches allowed is another pa­ rameter which controls the number of potential helices. The larger the value, the more the potential helices. In most ca~es we tested, allowing mismatches within 10% (e.g. 1 mismatch for 10 sequences) produces sa~isfactory results. However, the optimum value of the mtsmatches allowed will depend on how much variation exists be­ tween the sequences and how representative they are.

Weight values of the scoring function do not affect the number of potential helices at all, but has a direct and most influential effect on the prediction. Longer helices are more stable than shorter ones and exact-variant pairs are stronger evidence to base pairing than other types of pairs. Hence, the length and frequency of exact-variant pairs are assigned positive weights in the scoring func­ tion. Wobble pairs and inexact pairs are less stable than canonical pairs, and longer hairpin loops are less sta­ ble than shorter ones. Therefore, they have negative weights. While decision on the signs of the five param­ eters is straightforward, decision on the absolute values of them is not. Currently we do not have a general theo­ retic relation which explains how the relative sizes of the five values are correlated with respect to the prediction. The parameter values in Fig. 3 have been obtained exper­ imentally. We have successfully tested the program on different types ofRNAs (tRNAs, 5S rRNAs, 16S rRNAs, TAR, RRE, and cis-acting packaging sequences of HIV1), in many of which a widely used conventional program fails. A detailed comparison is described in [Han and Kim, 1993]. The execution time for predicting common foldings of 10 tRN As of length 77, 34 5S RN As of length 143, 10 16S rRNAs of length 1651, 10 TAR RNAs of length 57, 10 RRE RNAs of length 196, and 10 HIV-1 cis-packaging sequences of length 113 are 0.1, 0.6, 5.1, 0.1, 0.4, and 0.2 seconds, respectively.

5 Related Work

There currently exists no reliable and automatic method for modeling the folding process and for predicting folded structures. Existing approaches have limitations in their use either because they are computationally too expen­ sive to be practical, their predicted structures are not stable or they heavily rely on manual work [Brown, 1991; Noller and Woese, 1981; Zuker and Stiegler, 1981]. The algorithm by Sankof£[1985] for simultaneously aligning and folding M sequences of length N requires O(N 3M) time and O(N2M) space, and is not practical either in terms of time or space.

Like us, several researchers in AI have used approx­ imation in reasoning. Weld[1986] describes a simulator which uses aggregation as an abstraction technique for a continuous process. His approach cannot be applied to simulate the folding process because the process de­ pends on the constituent bases rather than on any reg­ ularity of the process, such as repeating cycles. Most research in qualitative kinematics has focused on pre­ dictin~ systems' behavior from the description of the sys­ tems lFaltings, 1987; Forbus et al., 1991; Gelsey, 1991; Joskowicz, 1989]. A general constraint common in those systems is that objects cannot occupy the same space. The RNA folding problem has a similar constraint in the sense that a base cannot pair with more than one base. However, the current technique of qualitative kinematics seems inapplicable to the RNA folding problem.

6 Future Work

The model presented in this paper has limitations and offers several directions for future research: (1) The pro­ cess of examining the effects of parameter changes of the model on the result is not automated. We plan to de­ velop methods, both theoretical and experimental, for automating this process. This will enable us to do more extensive perturbation study on the parameters affecting the behavior of the model. (2) We currently use a dis­ crete model both for the structure and for the structure formation process. For the structure, a discrete model is more appropriate than a continuous model because there are . no continuously varying parameters to consider in the model by its nature. The process of the structure formation is also modeled in a discrete model by taking "snapshots" of the structures at a series of intermediate stages. We might want to consider a continuous model for the process, which enables us to take into account the rate that a structure unwinds its helices to form a new helix as RNA grows. (3) The current model consid­ ers secondary interactions only and predicts a secondary structure. The model can be extended to incorporate tertiary interactions such as pseudo knots, which are base pairings between a loop and another part. Some of conflicting diagonals in the covariation matrix consti­ tute pseudo knots. ( 4) The current model can be easily applied to find disulfide bonds of protein structure as well. The model should be extensible to other general problems which involve matching or compatibility be­ tween objects. Any relation of compatibility or conflict between objects is represented in a covariation matrix, and a searching algorithm operates on the matrix to find all the possible solutions which satisfy the relation. The extended model is currently being implemented.

7 Conclusion

In this paper, we have described a new qualitative model for the folding problem of a biological system, and demonstrated that the model can be used both for pre­ dicting structures and for simulating the process of the structure formation over time. This model is a rough approximation of conventional models, but greater ac­ curacy in the prediction is obtained by using knowl­ edge from diverse sources instead of relying on a single source of knowledge, which is often inaccurate and in­ complete. Homologous RNA sequences are analyzed si­ multaneously for filtering potential helices and the result is incorporated with knowledge of free energy. A heuris­ tic algorithm with folding constraint propagator is used for searching for possible structures. The most plausible structures conserved in multiple homologous sequences are produced in the same amount of time and space as those of the conventional models for a single sequence. The model has been successfully tested on various types of RN As. It is known that there is no reliable or au­ tomatic way of predicting a common folding even if the related sequences are already aligned and that one way to obtain a common folding is to fold each sequence sep­ arately using an ordinary folding program and to search for common structures [Zuker, 1989]. Our model is the first one provably practical both in terms of time and space for finding folding structures of homologous RNAs.

Acknowledgments

We thank Mick Noordewier and Hong-Jin Kim for their contributions to this effort.