Nuclear Magnetic Resonance DP Dynamic Programming

Sampling methods are central to many algorithmic methods in structural RNA bioinformatics, where they are routinely used to identify important structural models, provide summarized pictures of the folding landscapes, or approximate quantities of interest at the thermodynamic equilibrium. In all of these examples, redundancy within sampled sets is uninformative and computationally wasteful, limiting the scope of application of existing methods. In this thesis, we introduce the concept of non-redundant sampling, and explore its applications and consequences in RNA bioinformatics.

We begin by formally introducing the concept of non-redundant sampling and demonstrate that any algorithm sampling in Boltzmann distribution can be modified into non-redundant variant. Its implementation relies on a specific data structure and a modification of the stochastic backtrack to return the set of unique structures, with the same complexity.

We then show a practical example by implementing the non-redundant principle into a combinatorial algorithm that samples locally optimal structures. We use this tool to study the RNA kinetics by modeling the folding landscapes generated from sets of locally optimal structures. These structures act as kinetic traps, influencing the outcome of the RNA kinetics, thus making their presence crucial. Empirical results show that the landscapes generated from the non-redundant samples are closer to the reality than those obtained by classic approaches.

We follow by addressing the problem of the efficient computation of the statistical estimates from non-redundant sampling sets. The absence of redundancy means that the naive estimator, obtained by averaging quantities observed in a sample, is erroneous. However we establish a non-trivial unbiased estimator specific to a set of unique Boltzmann distributed secondary structures. We show that the nonredundant sampling estimator performs better than the naive counterpart in most cases, specifically where most of the search space is covered by the sampling.

Finally, we introduce a sampling algorithm, along with its non-redundant counterpart, for secondary structures featuring simple-type pseudoknots. Pseudoknots are typically omitted due to complexity reasons, yet many of them have biological relevance. We begin by proposing a dynamic programming scheme that allows to enumerate all recursive pseudoknots consisting of two crossing helices, possibly containing unpaired bases. This scheme generalizes the one proposed by Reeders and Giegerich, chosen for its low time and space complexities.

III

We then explain how to adapt this decomposition into a statistical sampling algorithm for simple pseudoknots. We then present preliminary results, and discuss about extensions of the non-redundant principle in this context.

The work presented in this thesis not only opens the door towards kinetics analysis for longer RNA sequences, but also more detailed structural analysis of RNAs in general. Non-redundant sampling can be applied to analyze search spaces for combinatorial problems amenable to statistical sampling, including virtually any problem solved by dynamic programming. Non-redundant sampling principles are robust and typically easy to implement, as demonstrated by the inclusion of non-redundant sampling in recent versions of the popular Vienna package.

Résumé Substantiel

Des Acides RiboNucléiques, ou des ARNs, sont des biopolymères composés des nucléotides A, C, G et U. Ces séquences ont la propriété intéressante de se replier sur eux-mêmes, leur donnant la possibilité de former des structures différentes qui peuvent exercer des différentes fonctions. Les méthodes empiriques de leur déterminaison étant trop chères, lentes ou limitées, ces structures sont la cible d'étude du domaine de la bioinformatique structurale d'ARN.

Grâce au nombre des structures qu'une séquence moyenne d'ARN peut former, l'échantillonnage statistique est devenue une partie instrumentale à de nombreuses méthodes algorithmiques pour la bioinformatique structurale d'ARN, servant à identifier des structures secondaires importantes, des éléments clés d'espace de repliement suivant l'évolution de ces structures et estimation statistique des quantités d'intérêt de ces structures à l'équilibre thermodynamique. Des méthodes qui sont disponibles aujourd'hui souffrent d'effet de la redondance qui est en général non-informative et limite l'efficacité générale de ces algorithmes. Pour cette raison, dans cette thèse nous introduisons le principe de la génération nonredondante, permettant à éviter cet effet de manière optimale, et nous analysons ses applications dans le domaine de la bioinformatique structurale d'ARN.

Après avoir introduit le concept général de l'échantillonnage non-redondant, nous montrons que, à l'aide d'un formalisme que nous avons introduit, tout algorithme échantillonnant dans la distribution de Boltzmann peut être transformé en une version non-redondante. Ensuite nous expliquons le besoin de mémoriser des structures secondaires échantillonnées ainsi que le processus de leur sélection. Pour le faire, nous avons recherché une structure de données qui constitue une couche autonome à implémenter à l'algorithme qui est modifié en version nonredondante et que cette version est de même complexité que l'échantillonnage classique.

Pour une démonstration pratique, nous avons implémenté ce principe dans un algorithme de l'échantillonnage des structures secondaires localement optimales. Ces structures agissent comme des pièges cinétiques et impactent le résultat final du repliement, elles sont donc instrumentales lors de l'étude cinétique d'ARN. Nous alors utilisons cet outil pour étudier la cinétique d'ARN en modélisant l'espace de repliement des structures secondaires -espace représentant l'ensemble des structures secondaires et des trajectoires de repliement entre eux -et puis réalisant une intégration numérique sur ce modèle. Par la suite, des expériences mon-V trent que des modèles obtenus à partir des échantillons non-redondants présentent une cinétique plus rapide, et sont donc plus proches à la réalité, que ceux à partir des échantillons classiques.

Nous suivons par l'inclusion de l'échantillonnage non-redondant dans la version récente de la bibliothèque Vienna package populaire. Cette implémentation montre la facilité et la robustesse de l'implémentation du principe non-redondant dans des algorithmes existants dans le domaine de la bioinformatique structurale d'ARN. Ensuite nous montrons par des expériences que l'échantillonnage nonredondant de même algorithme est plus efficace que sa version de base et nous discutons sur l'importance d'optimisation de la structure de données sur la performance finale d'échantillonnage non-redondant.

Nous résoudrons par la suite le problème d'estimation statistique efficace à partir des échantillons non-redondants. L'estimation naïve, obtenu en moyennant des valeurs observées pour tous les échantillons, n'est pas possible dans ce cas car en absence de la fréquence, le résultat est biaisé. Par contre, nous sommes arrivés à établir un estimateur non-trivial et non-biaisé qui est spécifique aux échantillons non-redondants qui suivent la distribution de Boltzmann. Nous montrons sur des expériences concrètes sur des cas modèles que l'estimateur proposé par nous est plus efficace que l'estimateur naïf, et ceci notamment dans le cas ou la majorité d'espace de recherche a été échantillonné.

Finalement, nous introduisons un nouvel algorithme d'échantillonnage des structures secondaires contenant des pseudonoeuds du type H, le type le plus simple consistant de deux hélices s'entrecroisant. Ces éléments, même les plus simples sont typiquement omis pour des raisons d'efficacité, bien que beaucoup d'entre eux possèdent une grande importance biologique. Nous commençons par proposer un schéma de programmation dynamique qui permet à générer n'importe quel pseudonoeud composé de deux hélices comportant des bases nonappariées. Cet algorithme est une extension d'un algorithme existant de Reeders et Giegerich, choisi pour son efficacité temporelle est spatiale ainsi que la facilité relative d'implementation. Par la suite, nous expliquons comment construire le schèma de la programmation dynamique adaptable à un algorithme d'échantillonnage des structures secondaires avec des pseudonoeuds de type H. Après l'introduction des résultats préliminaires, nous expliquons comment il est possible à modifier cet algorithme pour inclure le principe de la génération nonredondante, et on discute sur des améliorations possibles de cet algorithme.

Des travaux présentés dans cette thèse ouvrent la porte vers l'analyse des séquences VI d'ARN plus longues et les études plus détaillées d'ARN en général. L'échantillonnage non-redondant peut être appliqué pour étudier l'espace de recherche de tous les problèmes combinatoires susceptibles à l'échantillonnage statistique, montrant sa grande versatilité, et il trouve des applications aussi aux des problèmes en dehors du domaine de la bioinformatique structurale d'ARN. VII VIII First, I would like to sincerely thank to my PhD Supervisor, Yann Ponty, for his continuous support during my PhD. His enthusiasm and dedication allowed me to learn a lot of new informations and acquire new interesting capabilities in scientific domain, but also new transversal skills. His patience and passion helped me to pull through the hard parts of my PhD, and to learn how calmly overcome them. Besides being the best supervisor I could ever wish for, Yann is also the close friend that also helped me to deal with the difficulties related to the life abroad, and taught me how to tackle such problems. His advices and encouragement helped me through my PhD and will stay helpful in the future to define and advance my professional path. I would also like to thank Hélène Touzet, my PhD Co-supervisor, for her support, giving me an inspiration to progress my work, and for her valuable help and insights on my thesis' manuscript. I also really appreciate her time to listen the rehearsals of my PhD defense along with my PhD director, and bringing the suggestions for improvements. She is also the person I know for the longest time, giving me the initial idea to follow this PhD, for which I am especially thankful. I thank the rapporteurs of my thesis' manuscript, Frédéric Cazals and Ivo Hofacker, for dedicating their time to read my thesis' manuscript and bring not only the detailed reports and necessary corrections, but also interesting ideas for future advancing of the project of my thesis. I would like to thank them, and also Samuela Pasquali, Adeline Pierrot and Aïda Ouangraoua, for accepting being members of my jury, reading again through my work, participating at my PhD defense and asking me the important questions about that bring useful insights to further develop the presented project.

I would like to thank the members of our laboratory, notably Mireille Régnier for welcoming me in the research group and the laboratory as well as explaining to me how it functions. I thank also many of past and current members of our research group. Namely Afaf, with whom I was discussing a lot about (not only) our PhD projects. Christelle, who helped me to remind myself of the statistic lessons I had in the past. Alice and Amélie, who guided me through the life of Ecole Polytechnique and also shared many stories they had when they were studying here. And Wei, Hua-Ting and Dimitri, for many, many talks we had. I thank also the members of the collaborating research groups in Vienna and Montreal. For the group in Vienna, I am would like to thank Ronny Lorenz, for his IX many ideas and many things I learned from him while we were working together. Maria, for discussions on our projects as well as interesting presentations during conferences. Andrea Tanzer, for her insights on my thesis' project and many interesting workshops she helped to organize and I had a chance to participate on. For the group in Montreal, I thank Jérôme Waldispühl and Mathieu Blanchette for their valuable help with my first presentation on international conference, and Roman, Carlos and Antoine for fun times during various conferences.

Finally, I would also like to thank my family, notably my beloved mum, Alena, and my grandma, Elena. Both fully supported me with my decision to stay abroad and work on my PhD. They have shown patience, compassion and understanding for most difficult moments of my PhD, and helped me to overcome these moments with ease. They gave me strength and resolve to follow my dreams. Without them, it would be impossible to me to finish my PhD and anything I have achieved thus far.

X

List of Figures

Preamble

Nucleic Acids are one of the most crucial macromolecule types for every living being. This class encompasses both Deoxyribonucleic Acids (DNA) and Ribonucleic Acids (RNA). The chemical composition of both acid types is quite similar; the main building brick in both cases is nucleotide composed from a phosphate group, a pentose sugar (2-Deoxyribose in DNA, Ribose in RNA) and a nitrogenous base [START_REF] Baker | Molecular biology of the gene[END_REF]. In the case of DNA, these bases are Adenine, Guanine, Cytosine and Thymine, famously abbreviated as A, G, C and T. While the first three can be found only in DNA, Thymine is exclusive to it and in RNA its role is substituted by Uracil (U) instead.

Unlike the chemical composition, the spatial structure structure of DNA and RNA is quite different. DNA is usually found in a double-stranded form arranged in a double helix, where the two strands are 'zipped' together via hydrogen bonds between the nitrogenous bases in an antiparallel manner (Figure 1.1). The orientation of each strand is determined by the sugars it presents -the 5' end of the strand directly branches from 5-th atom of 2-Deoxyribose. Since a given nitrogenous base has only limited number of pairing partners with which it can create sufficiently strong hydrogen bonds, the DNA adheres to base-pairing rules. RNA, on the other hand, is much shorter and can be found mostly in single stranded form, which is similarly oriented by Ribose (Figure 1.2). The base-pairing also applies to the RNA, though the pairs can be created within the same molecule or with other independent RNA, DNA or proteins [START_REF] Milman | The structure of a DNA-RNA hybrid[END_REF]. In the case of RNAs, the so-called non-canonical base pairs, not respecting the usual base pairing rules, 1 can also occur [START_REF] Das | Non-canonical base pairs and higher order structures in nucleic acids: crystal structure database analysis[END_REF]. Consequently, the combination of intramolecular interactions offers a variety of structures an RNA molecule can adopt. Each strand is a chain of nucleotides, connected by phosophodiester bond and composed of Phosphate (purple), 2-Deoxyribose (black) and a nitrogenous base. The base has four main variants -Adenine (A, green), Guanine (G, Yellow), Cytosine (C, Blue) and Thymine (T, brown). The two strands are held together by hydrogen bonds between nitrogenous bases following the base-pairing rules. The 5' end and 3' end are determined by the orientation of the 2-Deoxyribose.

The structure of DNA and RNA also projects into the functions they can perform. The stability of the double helix of DNA [START_REF] Tinoco | Nucleic acid structures, energetics, and dynamics[END_REF] allows it to serve as a cellular memory -it stores the genetic information within the cellular nucleus and that information is replicated and propagated on the successor during cell division.

In the case of RNA, people often think that its main function is the transfer of the information contained within the DNA to the molecular machinery that assembles the proteins -ribosomes [START_REF] Crick | Central dogma of molecular biology[END_REF]. However, besides these RNAs, also known as messenger RNAs (mRNA) there is a wide variety of its other types with equally wide array of functions other than information transfer; since these RNAs do not encode proteins, they are called non-coding RNAs (ncRNA). Their function can vary from structural, including but not limited to ribosomal RNA (rRNA) and transfer RNA (tRNA), to catalytic and regulatory functions such as the ncR-NAs involved in CRISPR-Cas9 system [START_REF] Rath | The CRISPR-Cas immune system: biology, mechanisms and applications[END_REF]; RFAM database currently stores an information about 2687 functional RNA families [START_REF] Kalvari | Rfam 13.0: shifting to a genome-centric resource for noncoding RNA families[END_REF]. The function of an RNA depends mainly on its structure, which is determined by the sequence and base-pairing rules, and consequently the structure of such sequence is often strongly conserved [START_REF] Williams | Phylogenetic analysis of tmRNA secondary structure[END_REF]. This, along with the functional versatility, motivates the studies aiming at modeling the structure of RNA. In this case, the sugar (black) is ribose and Uracil (U, orange) substitutes Thymine (T, brown).

The first methods of studying the secondary structures were empirical, which became standard over the years of their application. Such methods include nuclear magnetic resonance (NMR) [START_REF] Lincoln | RNA structure determination by NMR[END_REF], X-ray crystallography [START_REF] Westhof | Twenty years of RNA crystallography[END_REF] and cryo-electron microscopy [START_REF] Rees F Garmann | Visualizing the global secondary structure of a viral RNA genome with cryo-electron microscopy[END_REF]. Since such methods of finding an RNA structure are rather costly, complex and/or slow, with the advent of computational technologies it became tempting to exploit them to predict these structures. However, the task is not as simple as it may seem. First, when computing the structures the one with the lowest energy is considered the functional one. The reason for such assumption is that such structure is, under the conditions of thermodynamic equilibrium, where the system does not evolve anymore, considered the most stable, as it presents the global optimum. In reality however, this is not always the case. Instead, it might be a locally optimal structure where the RNA remained 'stuck' during the process of the creation of the structure, also known as RNA folding. These structures act as kinetics traps due to them being stabilized by sufficiently high energy barriers meaning high energy is needed to escape them [START_REF] Dotu | Computing folding pathways between RNA secondary structures[END_REF]. Another possibility is that the half-life of some RNAs is too short [START_REF] Steglich | Short RNA half-lives in the slow-growing marine cyanobacterium prochlorococcus[END_REF] and the molecule decays before it had enough time to reach the most stable structure. Second, some RNAs present more than one active structure, such as riboswitches, where the second structure is formed only in presence of certain ligand or during cotranscriptional folding when some of bases are not accessible [START_REF] Gong | Cotranscriptional folding and regulation mechanisms of riboswitches[END_REF]. This implies that merely studying the globally optimal RNA structure may not be sufficient to understand its relation with function. Rather, it is necessary to study also the so-called suboptimal structures -those with close but higher energy than the possible minimum, and also structural dynamics -simply stated the kinetics of RNA structures. Due to a number of structures that increases exponentially with the length of the RNA sequence [START_REF] Zuker | RNA secondary structures and their prediction[END_REF], exhaustive analysis in unrealistic, and it is necessary to select few representative structures that have high probability of being an active structure. This is the main motivation for the development of methods relying on statistical sampling.

While the decomposition algorithm that would become the basis and standard for the statistical sampling of secondary structures was discovered back in 1984, it took almost twenty years to establish the first stochastic sampling method. That was done by Y Ding and Charles E. Lawrence in 2003 [START_REF] Ding | A statistical sampling algorithm for RNA secondary structure prediction[END_REF]. These methods sample the structures by computing their probability according to the given distribution by decomposing these structures and sampling their local elements [START_REF] Li | Finding stable local optimal RNA secondary structures[END_REF]. Such sampling can be also done for a specific class of secondary structures, such as locally optimal ones.

Today, many sampling methods are available and can be used to retrieve the most interesting secondary structures according to given criteria. Unfortunately, most of them do not distinguish whether the structure they sample was generated before -they return it whether it was never sampled before, or it was already observed ten times. These repetitions are generally uninformative as shown by methods that partially remove them [START_REF] Kucharik | Basin Hopping Graph: a computational framework to characterize RNA folding landscapes[END_REF], and the time used to sample repeated structures might be used instead to sample the previously unseen examples that provide new information. The method sampling unique secondary structures would allow to cover the space of secondary structures much quicker. This can be used for faster and easier access to structures that have usually low probability of being sampled in the given distribution but might have biological significance, such as in the case of transitive structures [START_REF] Cupal | Density of states, metastable states, and saddle points: Exploring the energy landscape of an RNA molecule[END_REF], but also for the generation of a richer space of secondary structures to, for example, propose more accurate models of the RNA energy landscapes. For this reason, we asked whether there was a way to prevent such repeats, a method that is mainly, but not necessarily only, applicable in the domain of the structural RNA bioinformatics.

The main question is whether such approach can be performed in an efficient manner. Since the non-redundant sampling principle necessitates to perform supplementary computations and necessitates more information than the usual, redundant sampling approach [START_REF] Andrew | Non-redundant random generation algorithms for weighted context-free grammars[END_REF], it will obligatorily consume more time. Therefore, such method must be as optimal as possible so the gain from the nonredundancy outweighs the added cost, which must be reduced to minimum. It must be also easy to implement to all existing secondary structure sampling algorithms, meaning it should constitute a separate layer that is added to the algorithm to which the non redundant principle is implemented.

The work presented within this PhD answers all of the above questions and demonstrates the utility of such approach. Namely, we show how the non-redundant sampling can be used to access more precise prediction of the RNA energy landscapes, and to estimate more precisely the statistics related to their population.

The Plan of This Work

In this Section, we briefly remind the history of nucleic acids before introducing the basic concepts such as the levels of the structure. This is quickly followed by the Section 2, which elaborates on state-of-the art methods. After introducing notions and concepts necessary to understand this document as well as the notion of dynamic programming, we detail the algorithms employed in structural RNA bioinformatics that were crucial for our work. In this section, we also introduce a formalization of dynamic programming schemes. Such formalism helps us to generally demonstrate the properties of all algorithms that can be formalized and show the general compatibility with non-redundant sampling principle presented in the next chapter. Finally, we describe the state-of-the-art in RNA kinetics as well as the description of the modeling processes we employed.

We start Section 3 by describing the combinatorial algorithm that served as our starting point, called Algorithm of Saffarian [START_REF] Saffarian | RNA locally optimal secondary structures[END_REF], as well as justifying this choice. We explain its modification to compute the energies of structures it returns. We then pass on the principle of non-redundant sampling and how it can be implemented into the algorithm. The next part of the chapter describes the experiment of modeling the RNA kinetic landscape from the generated samples by this and competitor algorithms. We close this chapter by describing how the nonredundant sampling principle can be implemented into other existing algorithms while demonstrating it on the implementation into VIENNARNA library [START_REF] Lorenz | Vi-ennaRNA Package 2.0[END_REF], followed by the comparison of its performance with the basic version.

The Section 4 concentrates of statistical analysis of non-redundant sampling algorithm. The loss of the frequency due to non-redundancy means that the computation of an average of given property cannot be achieved naively by computing the weighted average. This problem was solved in collaboration with Christelle Rovetta by establishing a non-trivial, unbiased estimator. The details of it, as well as the experiment where we compare it with the naive variant, are described here.

In Chapter 5, we introduce the sampling algorithm for the secondary structures including H-type pseudoknots. This chapter explains the algorithm used in software pknotsRG [START_REF] Reeder | pknotsRG: RNA pseudoknot folding including near-optimal structures and sliding windows[END_REF], also extended and employed in pKiss [START_REF] Janssen | The RNA shapes studio[END_REF], to which our algorithm is an extension,as well as the reasons why we decided to extend this algorithm. We then present the preliminary results obtained as well as the future extensions, including using the non-redundant sampling principle.

In the Section 6, we summarize the contributions of non-redundant sampling principle and we discuss of their importance within and outside the field of structural RNA bioinformatics. This Section also serves as the conclusion, and includes a discussion of the future steps to advance the research presented here.

Nucleic Acids: A Look to the Past 1.3.1 The Discovery of Nucleic Acids

The discovery of nucleic acids, or more precisely of the DNA, dates back to 1869, when Friedrich Miescher discovered inside human leucocytes a substance that differed from proteins by its properties -he called it 'nuclein' [START_REF] Xie | DNA and the brains behind its discovery[END_REF][START_REF] Dahm | Friedrich Miescher and the discovery of DNA[END_REF]. The name was later changed by Richard Altmann to nucleic acid due to acidic properties the substance exhibited [START_REF] Dahm | Discovering DNA: Friedrich Miescher and the early years of nucleic acid research[END_REF]. Between 1885 and 1901, Albrecht Kossel and his students then identified five organic compounds present in nucleic acids; these later turned up to be nitrogenous bases. He also inferred that nucleic acids are most likely implicated in the synthesis of new tissues.

The differences between RNA and DNA were not apparent at the time of their discovery and at first they were called by their source of origin -DNA was known as the thymus nucleic acid [START_REF] Levene | On the structure of Thymus Nucleic Acid and on its possible bearing on the structure of plant nucleic acid[END_REF] and RNA as the yeast nucleic acid [START_REF] Levene | The structure of Yeast Nucleic Acid[END_REF], and if it was thought they were exclusive to animals and plants respectively. The key discovery was made in 1933 by Jean Brachet who has shown that the DNA was localized in chromosomes of every cell [START_REF] Brachet | Nucleic acids in the cell and the embryo[END_REF]. Likewise, he has shown that the RNA could be found in cytoplasm of every cell independently of its source.

Around the same time, Phoebus Levene studied the chemical composition of nucleic acids. His analysis revealed the chemical structure of each nucleotide, the nature of sugars as pentoses, the presence of phosphates and the fact that thymine and uracil are exclusive to DNA and RNA respectively [START_REF] Hargittai | The tetranucleotide hypothesis: a centennial[END_REF]. He also supposed that the DNA is composed from repeated tetranucleotide units interconnected by phosphate ester-sugar bonds. In 1944 Oswald Avery, Colin MacLeod and Maclyn McCarty performed a modified version of Griffith's experiment [START_REF] Griffith | The Significance of Pneumococcal Types[END_REF][START_REF] O T Avery | Studies on the chemical nature of the substance inducing transformation of pneumococcal types : Induction of transformation by a Desoxyribonucleic Acid fraction isolated from Pneumococcus Type III[END_REF]. In this version the mice were injected by rough (non-virulent) strain of Pneumococcus bacteria and the material of heat killed smooth (virulent) strain that was previously purified and that contained DNA. The results have shown that the DNA was the most likely responsible for bacterial transformation, and therefore it was involved in heredity. Consequently, it also disproved the tetranucleotide theory since it was unlikely that a mere chain of tetranucleotide repeats could contain all needed information.

The following important discovery would be made by Erwin Chargaff near the middle of 20 th century. He discovered two tendencies within the DNA that are today known as Chargaff's rules [START_REF] Elson | On the desoxyribonucleic acid content of sea urchin gametes[END_REF][START_REF] Chargraff | Composition of the desoxypentose nucleic acids of four genera of sea-urchin[END_REF]. The first discovery was that the content of Thymine and Adenine were equal, as well as that of Guanine and Cytosine. In other words %A = %T and %C = %G. This also definitely invalidated the tetranucleotide theory of Levene since not all nucleotides were found in equal proportions. The second was that these percentages were specific to species. From here, it was only one step towards the famous discovery made in 1953 by James Watson and Francis Crick, who using the X-ray crystallography images of DNA obtained by Rosalind Franklin and Maurice Wilkins created the two-strand double-helix DNA model [START_REF] Braun | How Rosalind Franklin discovered the helical structure of DNA: Experiments in diffraction[END_REF][START_REF] Watson | Molecular structure of nucleic acids. A structure for deoxyribose nucleic acid[END_REF]. This, along with first Chargaff's rule, has also shown that A pairs with T and C with G; consequently these pairs were called Watson-Crick's.

The boom around the DNA this discovery has caused then resulted into quick discovery of the genetic code in 1960s [START_REF] Nirenberg | The dependence of cell-free protein synthesis in e. coli upon naturally occurring or synthetic polyribonucleotides[END_REF][START_REF] Gardner | Synthetic polynucleotides and the amino acid code[END_REF][START_REF] Wahba | Synthetic polynucleotides and the amino acid code[END_REF] and the sequencing methods about dozen years later [START_REF] Sanger | A rapid method for determining sequences in DNA by primed synthesis with DNA polymerase[END_REF].

mRNA and ncRNA, and the discovery of the secondary structure

The discovery of DNA structure has shown that it is very unlikely that DNA itself would be directly transcribed to proteins. The experiments of Jean Brachet in forties of 20 th century suggested that the RNA is related to protein synthesis [START_REF] Brachet | Nucleic acids in the cell and the embryo[END_REF]. This led to Crick to postulate, in 1957, his Central Dogma of Molecular Biology [START_REF] Crick | Central dogma of molecular biology[END_REF] that suggested that RNA is the intermediate for passing the information from the DNA in nucleus to the center of protein synthesis in cytoplasm. This was later confirmed by the discovery of lac operon of E. Coli and experiments performed on them by Jacques Monod and Francois Jacob [START_REF] Monod | Teleonomic mechanisms in cellular metabolism, growth, and differentiation[END_REF]. The unstable RNA they found has become known as messenger RNA (mRNA) since then. This finding however opened a can of new questions like what the intermediate molecule that associates the genetic code to the correct amino acids is and of course the protein synthesis machinery itself. This led to the discovery of the first classes of the non-coding RNAs (ncRNA).

The ribosomes, the center of protein synthesis, are themselves constituted from a special subclass of ncRNAs -ribosomal RNA (rRNA). It was first observed in 1938 by Albert Claude [START_REF] Palade | Albert Claude and the beginnings of biological electron microscopy[END_REF] as a cytoplasmic masses he called microsomes, eventually identified as ribosomes by George Palade [START_REF] Palade | Studies on the endoplasmic reticulum. II. Simple dispositions in cells in situ[END_REF]. The transfer RNA (tRNA) -the link ensuring the correspondance between the genetic code and the aminoacids -was theoretically assumed by Crick. It was confirmed by Robert Holley et al when he reported a 77 nucleotides long sequence of yeast alanine tRNA along with suggestions of secondary structure for it in 1964 [START_REF] Holley | Structure of a Ribonucleic acid[END_REF]. The last point is, within the context of this thesis, particularly important because it demonstrates how ncRNAs, to which tRNA and rRNA belong, need specific structure to perform their function. It was then no surprise that a number of methods to identify these structures, first empirical like X-ray Crystallography and Nuclear Magnetic Resonance (NMR), later computational predictions, emerged. This thesis concentrates mainly on efficient computational prediction of RNA secondary structures and their kinetics.

A Structural Intermezzo -Levels of RNA structures

Before moving on towards the state-of-the-art in RNA energy computation and secondary structures prediction, it is necessary to precise what an RNA secondary structure is. That is the main objective of this subsection.

As already mentioned, RNAs consist of nucleotides (or ribonucleotides) A, G, C and U, which substitutes T from DNA. Each nucleotide consists of a nitrogenous base (from which the nucleotide gets its letter), Ribose and phosphate group. Nitrogenous bases can be also divided into Purines, denoted R, that have two heterocycles (Adenine and Guanine, see Figure 1.3), and Pyrimidines, with symbol Y, that have only one (others). A group of only nitrogenous base and Ribose is called ribonucleoside. Nucleotides are connected by phosphodiester bond between 3 rd and 5 th atom of Ribose 5-phosphate, from which the 5' end and 3' end orientation comes. Unlike DNA, RNA is mostly found in single-stranded form. Depending on the types of non-covalent interactions that are taken into an account, we distinguish four levels of structures of RNA. The first three are also exemplified on Figure 1.4

• RNA Primary Structure: No non-covalent bonds are considered. In other words, this is the sequence of nucleotides. An example of a primary structure is given on Figure 1.4A. This structural level is the most important for the mRNAs since they contain the information translated into proteins, but the higher-grade structures depend on it as well.

• RNA Secondary Structure: Only hydrogen bonds between nucleotides are taken into an account. Basically, the secondary structure is the set of base pairs that is present for given conformation of specific RNA. An example of such structure is given on Figure 1.4B. This level of structure is the main point of interest of this thesis.

The RNA secondary structure depends on a specific set of base pairs that is considered for given sequence. The principle of base pairs is hydrogen bond type of interaction between the nucleotides. Not every base pair can create these bonds. The usually considered base-pairs are those discovered by Watson and Crick in DNA -the pairs between A and U and C and G [START_REF] Watson | Molecular structure of nucleic acids. A structure for deoxyribose nucleic acid[END_REF]. The G -C pair is stronger [START_REF] Stofer | Free Energy calculations of Watson-Crick base pairing in aqueous solution[END_REF] due to, if we simplify enough, three hydrogen bonds instead of two that are formed between A and U (see Figure 2.1). Hence the stability of the RNA secondary structure is also considered by the number of G -C pairs, though by all means that is far from the only criterium. There are also other types of base pairs that are called 'wobble'. These pairs were predicted by F. Crick [START_REF] Crick | Codon-anticodon pairing: The wobble hypothesis[END_REF] due to the fact that while there are 64 codons, the number of corresponding tRNAs is lower [START_REF] Chan | GtRNAdb 2.0: an expanded database of transfer RNA genes identified in complete and draft genomes[END_REF], and consequently some of anticodon sites on some tRNA molecules have to pair with more than one nucleotide. Most of wobble base pairs include a special nitrogenous base, hypoxanthine, that is created by post-transcriptional modification of adenine in tRNA [START_REF] Tuorto | Genome recoding by tRNA modifications[END_REF]. Since that base is omitted within the work, we will consider, in addition to Watson-Crick base pairs, the only wobble base pair not including hypoxanthine, a base pair G -U (Figure 2.1, right).

In the RNA structures, there are also non-canonical base pairs that do not adhere to base pairing rules [START_REF] Das | Non-canonical base pairs and higher order structures in nucleic acids: crystal structure database analysis[END_REF]. We will not consider them here.

For the rest of this work, when talking about the secondary structures, we refer to the RNA secondary structures, unless noted otherwise.

• RNA Tertiary Structure: All intramolecular bonds are considered. That includes, besides hydrogen bonds between nucleotides, also hydrophobic effects and disulfide bonds. These forces give an RNA molecule a specific tridimensional structure. To have an idea of difference between secondary and tertiary structure, in the case of DNA, the latter is the typical doublehelix form, while its secondary structure would be merely a 'ladder'. An example of a tertiary structure is given on Figure 1.4C.

Like in the case of the secondary structures, there is an ongoing work on computational prediction and modeling of tridimensional RNA structures. Some are based on fundamental laws of physics, such as Molecular Dynamics, which describes the forces between different particles of RNA molecule by a set of parameters called the force fields [START_REF] Magnus | Computational modeling of rna 3d structures, with the aid of experimental restraints[END_REF]. Such methods are accurate but slow. Simplifications to such models can be made using experimental data, but these experimental data can be hard and expensive to obtain.

The RNA secondary structure is easier to predict than the tertiary structure due to the fact that the number of interactions to predict is lower. Nevertheless, it already provides a good basis for functional hypotheses, as exemplified by the regulation of tryptophan production on trp operon in E. Coli [START_REF] Du | trp RNA-binding attenuation protein-5' stem-loop RNA interaction is required for proper transcription attenuation control of the Bacillus subtilis trpEDCFBA operon[END_REF].

For this reason we decided to concentrate on the kinetics and prediction of RNA secondary structures. The higher level structures will not be developed further here and are only mentioned for the sake of completeness.

• RNA Quarternary Structure: In addition to the above, the intermolecular bonds are considered. This level of structure also requires the study of RNA-RNA, RNA-DNA and RNA-protein interaction. The secondary structure taken from RNA STRAND [START_REF] Andronescu | RNA STRAND: the RNA secondary structure and statistical analysis database[END_REF] and visualized by VARNA [START_REF] Darty | VARNA: Interactive drawing and editing of the RNA secondary structure[END_REF].

GCCCG GAUAG CUCAG UCGGU AGAGC AUCAG ACUUU UAAUC UGAGG GUCCA GGGUU CAAGU CCCUG UUCGG GCGCC A G C C C G G A U A G C U C A G U C G G U A G A G C A P C A G A C U U U U A A P C U G A G G G U C C A G G G U P C A A G U C C C U G U U C G G G C G C C A 1 76 A B C
The tertiary structure created by PYMOL [START_REF] Schrödinger | The PyMOL Molecular Graphics System[END_REF].

11

Chapter 2

RNA secondary structure prediction -What has been done

Notions related to the RNA secondary structures

Before we start dwelling deeper into the realm of RNA sequences and secondary structures, we must establish the basic notions that will be used in this work. We do not consider any other symbols outside the grammar Ω.

Definition 2.1.2 (Substring):

A substring of w w w from position i to position j

(j ≤ i ≤ j ≤ n) is denoted as w[i, j]. A nucleotide at position i, i ∈ [1, n] is marked as w[i], w[i] ∈ Ω.
RNAs treated here are strictly linear, with 5' end being on its left side. The nucleotides can form a different base pairs within the sequence.

Definition 2.1.3 (Base pair):

A base pair is a pair of nucleotides at position i and j, i < j. It is denoted by (i, j).

For given sequence, many different base pairs can be formed.

Definition 2.1.4 (Base pair set):

The set φ pair (w) φ pair (w) φ pair (w) is the set containing all possible base pairs (i, j) for a sequence w. Each base pair (i, j) must satisfy the following conditions:

• Each (i, j) must be Watson-Crick or wobble base pair, ie.

(i, j) ∈ {{A, U}, {C, G}, {G, U}}.

• w[i] and w[j] must be separated by at least θ base pairs:

j -i -1 ≥ θ.
The minimum base pair length comes from the steric strain that forbids the base pairing of nucleotides way too close one to another.

Each sequence w can form a certain subset of all possible base pairs φ pair (w) at any given moment. This applies to its substrings w[i, j] as well. From this point on, we refer an to RNA secondary structure as a secondary structure or just structure unless noted otherwise.

Throughout this work, we may refer to the secondary structures smaller than n and with undefined limits.

Definition 2.1.6 (Secondary structure with unspecified length): A secondary structure s is a structure with unspecified, resp. implicitly specified w[i], w[j] and n.

Within the same S(i, j), each base pair has to follow, in addition to the requirements from the Definition 2.1, the following conditions:

• One nucleotide can participate at most in one base pair at the same time.

This means that if (i, j) ∈ S, ∀k = j, (i, k) / ∈ S, (j, k) / ∈ S.

• At first, we consider for any base pairs (i, j) ∈ S, (k, l) ∈ S with i < k that either i < k < l < j or i < j < k < l, ie. one base pair cannot cross another. At later point, in Section 5, we will also investigate some structures that include crossing base pairs.

If S = ∅ (ie. all of bases within the sequence are unpaired), the RNA sequence is called unfolded.

Definition 2.1.7 (Space of secondary structures):

Given an RNA sequence w, we define by S, the space of secondary structures, a set of all possible secondary structures of any length. The space S n is a space of n-sized secondary structures observed for w w w. Each is composed of base pairs from φ pair (w). Similarly, S(i, j) the space of all secondary structures S(i, j) for w[i, j].

Representation of RNA secondary structures

A secondary structure can be shown in various ways. Each of them has its advantages and weak points. The representations used throughout this work are listed here along with an example on Figure 2.2.

• Radial representation (2.2A): The paired segments are represented by a rectangular ladder-like structure while the unpaired segments are organized in circles from which the paired regions branch. In the case of a linear RNA, its unpaired base can be represented linearly or in a circle. This visualization clearly depicts different local substructures (see Section 2.4.2), facilitating to understand its function. On the other hand, it badly handles the cases where base pairs cross. It is good to be used when one needs to explain the function of a secondary structure.

• Circular representation (2.2B): The (linear) RNA is represented in a circular pattern, base pair being represented by a line or an arc. It has no problem to depict crossing base pairs but it is much harder to interpret. This type of secondary structure will not be employed in this work.

• Linear representation (2.2C): The RNA is shown as a line, and base pairs are depicted by arcs connecting the concerned nucleotides. This representation can be used to compare two different structures each drawn from other side of line. It is also an easy way to depict and explain the recursions of secondary structure prediction algorithms, which is main use of this visualization here. However, it might be a bit harder to interpret, especially for long-range interactions.

• Dot-bracket notation (2.2D): This is the textual representation of secondary structure with specific symbols. A dot '.' indicates an unpaired base, and a base pair is represented by brackets '(' and ')'. The correspondence of opening and closing brackets works in the same manner as in math; six opened brackets must be closed by six closing brackets and the last bracket opened it the first one that is closed. In the case of crossing base pairs, one of the crossing base pairs can be denoted by square bracket possibly followed by a letter, ie. '['-']' or '[A'-']A'. This format is hard to interpret by human, though there is a visualization software that can depict the related secondary structure, such as VARNA [START_REF] Darty | VARNA: Interactive drawing and editing of the RNA secondary structure[END_REF]. Its main purpose is to be machine-readable, which is why it became a standard input and output format for RNA secondary structure prediction software. The Guanosine-Uridine (right) base pair is also called wobble and is less frequent due to slightly lower energy.

Combinatorics of RNA secondary structures

First methods used to find the RNA secondary structure and structures in general were empiric. Even the discovery of double helix of DNA resulted from an X-ray The different elements are an external loop (cyan), hairpin (red), stem/helix (blue), an internal loop (grey), a bulge (white) and a multibranched loop (yellow). Representations A-C were created using VaRNA [START_REF] Darty | VARNA: Interactive drawing and editing of the RNA secondary structure[END_REF].

A G C G C G G A A A A C G C G A A A A A A A C G C G A A A A C G C G A A A C G C G A A A A A A A A C G C G A C G C G A C G C G A A C G C G A A A A A A A A A C G C G A A C G C G A G C G C G A U A 1
A G C G C G G A A A A C G C G A A A A A A A C G C G A A A A C G C G A A A C G C G A A A A A A A A C G C G A C G C G A C G C G A A C G C G A A A A A A A A A C G C G A A C G C G A G C G C G A U A
A G C G C G G A A A A C G C G A A A A A A A C G C G A A A A C G C G A A A C G C G A A A A A A A A C G C G A C G C G A C G C G A A C G C G A A A A A A A A A C G C G A A C G C G A G C G C G A U A
crystallography picture [START_REF] Braun | How Rosalind Franklin discovered the helical structure of DNA: Experiments in diffraction[END_REF]. However the crystallography is impossible if the analyzed substance does not crystallize. The NMR does not have similar problems, and unlike for proteins, the length of the molecule is not a limit either due to the RNA being shorter molecules. However, NMR and other experimental methods to determine the structure of RNA such as the X-ray crystallography and cryogenic electron microscopy have complicated protocols that require some time to prepare and they are expensive even today [START_REF] Bottaro | Towards de novo RNA 3D structure prediction[END_REF]. Therefore with the advent of computational technology the attention turned to the computational prediction of the molecular structure, where the possibilities to make it cheap and more efficient are much larger.

The number of secondary structures

To predict a secondary structure, it is necessary to have a mathematical or combinatorial description of properties that such a structure has. This property has to be possible to be evaluated by a score, computed by appropriate scoring function S.

Definition 2.3.1 (Scoring function):

We call a scoring function a function that to each S ∈ S, associates a certain score:

S : S → R
In the case of RNA, S usually computes a free energy. The lower energy of the structure means the better stability of the final secondary structure and hypothetically this one is the most likely structure to be the functional one.

The first approaches of energy computation were simply based on a number of base pairs within the considered secondary structures, without any distinction between stronger G-C, weaker A-U and wobble G-U pairs. The decomposition method -necessary to establish an evaluation function -was introduced by Temple F. Smith and Michael S. Waterman [START_REF] Michael | RNA secondary structure: A complete mathematical analysis[END_REF][START_REF] Michael S Waterman | Secondary structure of single-stranded nucleic acids[END_REF] to compute the number of possible secondary structures. In principle, each base is, under the restriction defined in Section 2.4.2, either unpaired, or paired to some other base. With this in mind, if we suppose a substring w[i, j], then j is either unpaired, or paired to some k, i ≤ k < j -θ, θ implying the minimum base pair length. In the first case, the number of the secondary structures is equal to those of the substring w[i, j -1]. In the second case the pair (k, j) generates two new intervals

w[i, k -1]
and w[k+1, j-1] (see Figure 2.3). The same principle is then reapplied recursively to each new interval until j -i -1 ≤ θ.

Let N d S(i,j) be the number of secondary structures of w[i, j] having exactly d pairs. The recursive relation given by Smith and Waterman algorithm is:

N d+1 S(i,j) = N d+1 S(i,j-1) + j-θ-1 k=i d u=0 N u S(i,k-1) × N d-u S(k+1,j-1) (2.1)
Here the first component of the right-hand side part of this equality is the case where i is unpaired, and the second component, with the double sum, is the case where j is paired with some k. Note that for w[i, j] | j -i -1 ≤ θ no base pair is possible and every nucleotide must stay unpaired. Therefore,

∀{i, j | j -i -1 ≤ θ}, N 0 S(i,j) = 1 and ∀{i, j, d, j -i -1 ≤ θ, d > 0}, N d S(i,j) = 0.
Let N S(i,j) = card(S n (i, j)) be the number of secondary structures of w[i, j], or in other words, N S(i,j) = n d=0 N d S(i,j) , n being an upper limit of d. In a similar manner, by generalizing the decomposition from Figure 2.3, we get:

N S(i,j) = N S(i,j-1) + j-θ-1 k=i N S(i,k-1) × N S(k+1,j-1) (2.2) Here ∀{i, j | j -i -1 ≤ θ}, N S(i,j) = 1.
The computation of the number of the secondary structures is then simple. First we pre-compute all values of N S(i,j) 1 ≤ i < j ≤ n in a manner the the shorter sequences are precomputed first and store them in a two-dimensional matrix. The final value is stored within N S (1,n) . The time complexity of the algorithm is O(n 3 ) since n 2 entries must be treated and for each n options exist at most. The space complexity is O(n 2 ) due to storing n 2 entries. Both cases produce shorter substrings on which this process can be recursively applied.

i j i j-1 j i k-1 k k+1 j-1 j

About Dynamic Programming

Enumerating all possible secondary structures by simple recursive function is not the most efficient way to handle this problem. Using the algorithm of Smith and Waterman, Michael Zuker and David Sankoff shown that for the sequence with uniformly distributed nucleotides the asymptotic limit for the number of secondary structures of RNA of length n is around 1.867 n [START_REF] Zuker | RNA secondary structures and their prediction[END_REF] -it raises exponentially. Therefore, trying all possible combinations would result in an exponential complexity. This is why the approaches using Dynamic Programming (DP) are extremely useful. DP scheme decomposes a central problem, here the secondary structure of w, into smaller problems, here the structures of w[i, j], computing those and using the results to solve the main problem. Each calculation is done only once and the result stored in DP matrices whose size depends on the number of states that have to be saved, like done by counting function from Equation 2.1. This makes the computation much faster than a brute-force enumeration of all possible solutions, by avoiding the repeated computation of solution to a given problem.

The DP algorithms have to respect one main principle, called Bellman's Principle of Optimality [START_REF] Bellman | The theory of dynamic programming[END_REF]. To quote Bellman:

"An optimal policy has the property that whatever the initial state and initial decision are, the remaining decisions must constitute an optimal policy with regard to the state resulting from the first decision."

Translated, every intermediate solution, that leads to the main optimal solution, must be also optimal. These are optimal solutions to given subproblems. Note that the DP schemes might be also applied in cases where the algorithm would give an optimal solution from non-optimal ones for given subproblems or inversely, where optimal solution to subproblem does not lead to an optimal global solution, but that reduces them to a mere heuristic method.

All of the approaches used for structure prediction and mentioned here respect Bellman's Principle of Optimality unless noted otherwise. The core of most of them is to decompose w into w[i, j] | 1 ≤ i < j ≤ n and solve the subproblems for sequences of increasing length l = j -i + 1. The result for each subproblem is computed and memorized. Consequently, the results for substrings of length 1, 2...l -1 are used to compute that for sequences of length l. The main problem is then solved for the sequence w. The same principle can be employed to compute suboptimal structures. The equivalent definition can be established when maximizing S.

Nussinov Algorithm

The first algorithm to make use of the decomposition of Smith and Waterman, as well as the DP, to predict a secondary structure was proposed by Ruth Nussinov et al [START_REF] Nussinov | Algorithms for loop matchings[END_REF]. The function S is in this case the number of the base pairs.

The first version of the algorithm was not unambiguous, meaning that a single structure could be decomposed in multiple ways. Here we focus on an unambiguous version which was also the base for the secondary structure counting function of Waterman [START_REF] Michael | RNA secondary structure: A complete mathematical analysis[END_REF][START_REF] Michael S Waterman | Secondary structure of single-stranded nucleic acids[END_REF].

Here each base pair has the same weight; G -C, A -U and G -U score as 1. The final score N(1, n) gives the maximum number of the base pairs possible for w.

Let N i,j be this number for any w[i, j]|1 ≤ i < j ≤ n, then we can formulate the algorithm as:

N i,j = max        N i,j-1 if j is not paired max i≤k≤j-θ-1 (i,j)∈φ pair (w) (N i,k-1 + N k+1,j-1 + 1) if j is paired to k (2.
3)

The score raises by 1 whenever k and j are paired, otherwise the score remains unchanged. It can be generalized by replacing the bonus 1 in the case of (k, j) ∈ S(i, j) by a bonus B k,j that varies depending on a base pair type:

N i,j = max        N i,j-1 if j is not paired max i≤k≤j-θ-1 (i,j)∈φ pair (w) (N i,k-1 + N k+1,j-1 + B k,j ) if j is paired to k (2.4)
N i,j = 0 if j ≤ i + θ since the substring w[i, j] is too short for any base pair to be created there, providing an initialization condition. From there, we can compute N i,j for progressively longer substrings. The values are stored in a matrix N and reused as needed. The final result is stored in N 1,n .

To obtain the structure with the maximum base pairs, we can either directly store the base pairs as the maximum values are selected for different points, or perform a stochastic backtrack after computing N 1,n . The backtracking is done by starting at N 1,n , finding the combination that gives us the stored value -in this case it is either

N 1,n-1 or sum of N 1,k-1 and N k,n for 1 ≤ k ≤ n -θ -1. If (k, j) is present,
it is added to the structure, and the whole process is repeated until the initial conditions are reached. The secondary structure generated this way is the one with maximum base pairs.

The complexity of this algorithm is O(n 3 ) in time and O(n 2 ) in space [START_REF] Nussinov | Algorithms for loop matchings[END_REF]. The explanation is equivalent to the secondary structure counting function of Waterman.

Further extensions added no-lone-base-pair as a constraint -no base pair can appear alone. This is because the stacks of a base pair have stronger stabilizing effects on the structure than a single base pair due to the stacking effect of the cycles of nitrogenous bases [START_REF] Petersheim | Base-stacking and base-pairing contributions to helix stability: thermodynamics of double-helix formation with CCGG, CCGGp, CCGGAp, ACCGGp, CCGGUp, and ACCGGUp[END_REF]. Despite this, the prediction results of this model were still imprecise, and more complex approaches were necessary.

Enter the Thermodynamics

The approaches in question targeted to compute the exact free energy of the entire secondary structure. The free energy can be described, in the context of secondary structures, as the amount of the energy you would get/supply to the unfolded RNA sequence for it to reach the given secondary structure. Its unit is kcal.mol -1 , the reference being an unfolded state. The most interesting structure S for given w is the one with the lowest possible energy, called Minimum Free Energy structure (MFE).

To compute MFE, S must be decomposed into a set of elements for which the free energy can be computed without too much problems. Their free energy contribution is then summed up, netting the final result. This also implies that the elements must be independent one from another, or depend on a context which can be concisely described. The decomposition and the resulting elements depend on the energy model which will be applied to compute the energies.

About energy models

The energy model is a model that computes the free energy of a set of specific elements. These may vary from simple values such as the energy of G -C bond to much more complex relations that apply only for very specific cases, such as loops that have exactly four unpaired nucleotides, dangling ends -interactions between the base pair and its unpaired neighbors -and so on. Throughout this work, we will mainly talk about and employ the Turner Energy Model [START_REF] Turner | RNA structure prediction[END_REF][START_REF] Douglas | NNDB: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure[END_REF]. This energy model allows to compute the energy of the elements resulting from the decomposition of a secondary structure into loops, as detailed in Section 2.4.2.

The Turner Energy Model consists of values that are tabulated and were directly deduced from experimental results for smaller elements, such as the stabilities of terminal wobble G -U pairs in consecutive base pairs [START_REF] Jonathan L Chen | Testing the nearest neighbor model for canonical RNA base pairs: revision of GU parameters[END_REF]. The mathematical relations for bigger substructures with more nucleotides were extrapolated from the empiric results such as melting experiments of the RNA [START_REF] Susan | Optical melting measurements of nucleic acid thermodynamics[END_REF]. The energy of each of substructure, when calculated by this model, depends on and only on the properties of the substructure itself and of those of the adjacent substructures such as neighboring base pairs. For this reason, this model is also called Nearest Neighbor Turner Energy Model. This is an important property because it suggests the independence between the distant elements resulting from the decomposition and therefore the energy contribution of each of them can be simply summed up. This is crucial for the design of a dynamic programming algorithm.

All successive versions of this model, including its latest installment contributed in 2014, were implemented into VIENNARNA [START_REF] Lorenz | Vi-ennaRNA Package 2.0[END_REF], which is the library we use for the application of this energy model.

A decomposition of RNA secondary structure

There are multiple ways to decompose a secondary structure. The most simple way was presented by Nussinov algorithm, where the only criterium is whether the base is paired or not. For the Turner Energy Model to apply, the structures are decomposed into an elements called loops. These can be sorted into groups depending on their properties, and for which the energy is computed in the same manner. This way the DP schemes can be broken down into a number of cases depending on the element that is treated, similarly to distinguishing paired and unpaired nucleotides in the scheme used in Nussinov Algorithm. Due to the independence condition, the free energy E of a secondary structure S is given by:

E = l∈S E l (2.5)
where l stands for loop and E l is computed by Turner energy model for l that depends on it and potentially on neighboring substructures. The different loops are an external loop (cyan), hairpin (red), stem/helix (blue), an internal loop (grey), a bulge (white) and a multibranched loop (yellow). The secondary structure was created using VaRNA [START_REF] Darty | VARNA: Interactive drawing and editing of the RNA secondary structure[END_REF].

A G C G C G C A A A A C G C G A A A A A A A C G C G A A A A C G C G A A A C G C G A A A A A A A A C G C G A C G C G A C G C G A A C G C G A A A A A A A A A C G C G A A C G C G A G C G C G A U A
With exceptions, all loops are delimited by two types of base pairs. A single loop can contain different branchings. The base pair of a branching that leads to the 5' end and 3' end of the sequence is called a closing base pair, here denoted by (i, j), since it encloses the said loop. The pairs enclosed by loop themselves are called opening base pairs due to them opening other loops. These are not limited to one and they open the following substructures. Here, they are denoted (k x , l x ).

The opening base pair for a given loop is also a closing base pair for the following one. The decomposition is illustrated on Figure 2.4 and shows an example of each using different colors.

Here we list all possible loops that can appear in secondary structure S. Each loop has its specific energy to be computed by the Turner Energy Model (See Section 2.1 for definitions).

• Stems have stabilizing effect on the secondary structure, therefore a single base pair is a rather rare occurrence. Its energy is denoted by E St (i, j).

• An internal loop. A set of two strands of an unpaired nucleotides (i, j) and (k, l). The two unpaired strands must have the same length, or j -l = k -i (grey on Figure 2.4).

• A bulge. Similar to the internal loop, except the unpaired stretches are not of the same length α and β (white color on Figure 2.4). The energy by the Turner Energy Model is computed in the same manner for both internal loop and bulge, so for both cases we note E IL (i, j, k, l). In some cases a stem is also considered as a specific case of an internal loop with α = β = 0. In that case we note E ILG (i, j, k, l) such that

E ILG (i, j, k, l) =    E St (i, j) if k = i + 1 and l = j -1 E IL (i, j, k, l) otherwise
• A multibranch loop. Multiple stretches of unpaired nucleotides delimited by single closing base pair (i, j) and at least two opening base pairs (k x , l x ), x = 1..a, a ≥ 2 (yellow on Figure 2.4).

The energy computation of multibranch loops of the secondary structure is the most delicate part that would result into an exponential complexity of the associated DP scheme. An energy of a multibranch loop, here denoted E M (i, j), depends on many factors such as the types of opening and closing base pairs and the number of branchings along with other dependencies. However, it can be reasonably approximated by a simplified Jacobson-Stockmayer expression [START_REF] Hower | Parametric analysis of RNA branching configurations[END_REF][START_REF] Jost | Prediction of RNA multiloop and pseudoknot conformations from a lattice-based, coarse-grain tertiary structure model[END_REF]:

E M (i, j) = a + bN branch + cN unp (2.6)
where a is the penalty for the creation of the multibranch loop, b for each unpaired base it contains, c the penalty for each branch within it, and N unp and N branch the number of unpaired bases and branches within the multibranch loop respectively. The values of a, b, c are supposed to be constant. This way it is possible to decompose a multibranch loop into separate helices and compute E M (i, j) by summing up all contributions/penalties. For this reason most of MFE prediction algorithms employing Turner Energy Model use such linear approximation.

We also mention the type of an element that is not treated here due to the restrictions on structures but is addressed later in this work -a pseudoknot.

Definition 2.4.2 (Pseudoknot):

A pseudoknot is a structure that contains at least two pairs (i, j) and (k, l) such that i < k < j < l.

In other words a pseudoknot consists of crossing base pairs. There are different types of pseudoknots of increasing complexity [START_REF] Christian M Reidys | Topology and prediction of RNA pseudoknots[END_REF]. Due to the fact the introduction of pseudoknots into a secondary structures vastly increases their number and searching for the pseudoknots with the vastest definition is proved to be an NP-complete problem [START_REF] Akutsu | Dynamic programming algorithms for RNA secondary structure prediction with pseudoknots[END_REF][START_REF] Rune | Pseudoknots in RNA secondary structures[END_REF], they are omitted in the first part of this thesis. Later, we focus on an algorithm that treats an H-type class of pseudoknots, a problem that can be reduced to a polynomial complexity.

Zuker Algorithm

The first algorithm that made the use of this decomposition alongside the Turner Energy Model was introduced by Michael Zuker et al [START_REF] Zuker | Optimal computer folding of large RNA sequences using thermodynamics and auxiliary information[END_REF][START_REF] Zuker | RNA secondary structures and their prediction[END_REF][START_REF] Rune B Lyngsø | An improved algorithm for RNA secondary structure prediction[END_REF]. Unlike Nussinov algorithm, it looks for minimum free energy secondary structure. The decomposition of a secondary structure S in loops necessitates the scoring function S being decomposed in three values for every w[i, j]:

• W i,j -the minimum possible energy for S(i, j)

• V i,j -the minimum possible energy for S(i, j) | (i, j) ∈ S(i, j)

• WM i,j -w[i, j] is contained within a multibranch loop and contains at least one helix

The algorithm itself is given by:

W i,j = min    W i,j-1 min i≤k≤j-θ-1 (W i,k-1 + V k,j ) V i,j = min                    E H (i, j) E S (i, j) + V i+1,j-1 min i<k<l<j k≤l-θ-1 (E IL (i, k, l, j) + V k,l ) min i<k<j (WM i,k-1 + WM k,j + a) WM i,j = min                WM i,j-1 + c WM i+1,j + c V i,j + b min i<k<j (WM i,k + WM k+1,j ) . (2.7)
where a, b, c are constants and E H , E S , E IL are loop energies as defined in Section 2.4.2. For W i,j , w[j] is either unpaired or paired to some w[k] | i ≤ k ≤ j-θ-1. For V i,j , either every w[k] | i < k < j is unpaired and (i, j) forms hairpin, (i+1, j-1) ∈ S(i, j) forms a stack, (k, l) ∈ S(i, j) with i < j < k < l creates a bulge or internal loop or w[i, k-1] and w[k, j] contain each at least one base pair and therefore form a multibranch loop. VM i,j might have either w[i] or w[j] unpaired (i, j) ∈ S(i, j) or it can be split into two another helices. For VM i,j , we add the penalty c for each unpaired base found, and b every new base pair. The entire decomposition is illustrated on Figure 2.5.

The initialization condition is is ∀{i, j | j ≤ i + θ}, W i,j = V i,j = WM i,j = 0 due to w[i, j] being too short. Each step solves and stores the results of W i,j , V i,j and VM i,j in three separate DP matrices by increasing length. At the end, the MFE is given by W 1,n , the precise secondary structure being obtainable either by memorizing the pairs that contributed to the MFE or by traceback from W 1,n .

W i,j

i j i j-1 j i k-1 k j V i,j i j i j i i+1 j-1 j i k l j i i+1 k k+1 j-1 j W M i,j i j i j-1 j i i+1 j i j i k k+1 j Figure 2.5:
The principle of Zuker Algorithm. W i,j , V i,j and WM i,j represent the minimum energy of a structure S a substring w[i, j] can form, the minimum energy of S(i, j) can form if (i, j) ∈ S(i, j), and a minimum energy section of a multibranch loop with at least one helix respectively. W i,j : for a given substring w[i, j], j can be either unpaired or paired to some k such that i ≤ k ≤ j -θ -1. V i,j : if (i, j) exists, it can enclose only unpaired bases, a pair (i + 1, j -1), a single or multiple base pairs. VM i,j : The first of last base can be unpaired, (i, j) ∈ S(i, j) or it can contain more helices, in which case the substring is split in two by k.

One of major tools used for secondary structure prediction and structural analysis, VIENNARNA library [START_REF] Lorenz | Vi-ennaRNA Package 2.0[END_REF], uses a variant of Zuker algorithm. The modification touches notably multibranch loops, where a new scoring function M 1 i,j is introduced and WM i,j is replaced by M i,j . While the first M 1 i,j contains exactly one helix of multibranch loop within w[i, j], M i,j can contain more of them (but at least one). The algorithm itself, also depicted on Figure 2.6 is follows below.

F i,j = min    F i,j-1 min i≤k≤j-θ-1 (F i,k-1 + C k,j ) C i,j = min              E H (i, j) min i<k<l<j k≤l-θ-1 (E ILG (i, j, k, l) + C k,l ) min i<k<j M i+1,k-1 + M 1 k,j-1 + a + b M i,j = min      min i<k<j c × (k -i) + M 1 k,j min i<k<j M i,k-1 + M 1 k,j M 1 i,j = min M 1 i,j-1 + c C i,j + b (2.8)
Here the multibranch loop get split into the last helix M 1 k,j and the rest M i,k-1 .

M i,k-1 can contain either multiple helices or one helix preceded by a stretch of unpaired nucleotides. M 1 k,j might contain unpaired bases at its w[j] end. F i,j and C i,j are equivalent to W i,j , respectively V i,j from the base version. The initialization remains unchanged, ie. ∀i, j, j ≤ i + θ, F i,j = C i,j = M i,j = M 1 i,j = 0. F 1,n gives the free energy of the MFE for w.

The complexity of both algorithms is O(n 4 ) in time and O(n 2 ) in space. The limiting step is the evaluation done for the internal loops/bulges, which has four different parameters i, j, k, l, hence the O(n 4 ) complexity. This complexity can be reduced to O(n 3 ) if the maximum length of its unpaired stretches α and β is limited [START_REF] Mccaskill | The equilibrium partition function and base pair binding probabilities for RNA secondary structure[END_REF], though some modifications succeeded to reduce it without said restriction using some further assumptions on the energy model within internal loops [START_REF] Rune B Lyngsø | Internal loops in RNA secondary structure prediction[END_REF].

In the rest of this work, we will reference this version of the algorithm as VZu (VIENNARNA-Zuker) algorithm. and C i,j are functionally identical to W i,j , respectively V i,j from the original version of Turner's Algorithm. M i,j and M 1 i,j score the section of a multibranched loop within w[i, j] containing at least one helix and exactly one helix respectively. M i,j is split to a section that contains exactly one helix possibly preceeded by unpaired bases and another that is either completely unpaired or contains at least one another helix. M 1 i,j contains one helix potentially followed by unpaired bases.

F i,j i j i j-1 j i k-1 k j C i,j i j i j i k l j i i+1 k-1 k j-1 j M i,j i j i k-1 k j i k-1 k j M 1 i,j i j i j-1 j i j

From RNA thermodynamics to kinetics

MFE computations allow us to assess which structure is the most stable. However, it may not necessarily correspond to the empirically observed functional structure. This may result from the approximations made in Turner Energy Model to simplify the computation, such as approximating the free energy of multibranch loops in Equation 2.6. This may cause the MFE to be a different structure in reality, though usually with the similar energy.

Another possibility is that before MFE can form, the RNA sequence stays trapped in secondary structure that has lower energy that the neighboring structures. Definition 2.5.1 (Neighboring structure): We call a neighboring secondary structure of S a structure that differs from S by one base pair, including base pair shifts or for breaking and recreating a base pair with at least one common nucleotide. The space of neighbors of S is denoted as V S .

The definition of locally optimal structure follows. Definition 2.5.2 (Locally optimal structure): We call a locally optimal structure a structure S for which

S(S) < min

Sx∈V S (S(S x )).

These structures can act as potential kinetic traps, with significant impact on the process of the formation of secondary structure.

Finally, one sequence can have more functional structures, as in the example of riboswitches [START_REF] Jacob T Polaski | Cobalamin riboswitches exhibit a broad range of ability to discriminate between methylcobalamin and adenosylcobalamin[END_REF], so finding a single MFE will not find the so-called metastable structure. Sometimes the metastable structure results from co-transcriptional folding, ie. the structure starts to fold before the transcription process is completed, meaning some bases are not accessible at the beginning of the folding [START_REF] Gong | Cotranscriptional folding and regulation mechanisms of riboswitches[END_REF].

For these reasons the scientific community became interested in suboptimal secondary structures. The main interest of this work is to propose a method that allows us to retrieve the suboptimal structures in an efficient manner. The utility of suboptimal structures is they can be used to create a representation that is very useful for the study of the kinetics of secondary structures, an RNA folding landscape.

RNA folding landscape

Definition 2.5.3 (RNA Folding landscape): An RNA folding landscape is the model of the space S n of secondary structures S n for given sequence w, visualizing their energies and eventually the transitions between them.

It can be seen as a visualization of the system constituted of sequences w folded in every secondary structure S ∈ S n . Some of these structures can change one to another, which is represented by a transition within the folding landscape. It can also model only a subset S s sub ⊂ S n . Due to the complexity and multidimensionality of RNA folding landscape, it is difficult to depict the entirety of the landscape at once for longer sequences.

For this reason, most representations therefore depict only some section of the folding landscape depending on the used reference. For example the Figure 2.7 shows an example of two energy landscapes for Yeast Phenylalanine tRNA (Figure 2.7 left) [START_REF] Westhof | Restrained refinement of the monoclinic form of yeast phenylalanine transfer RNA. Temperature factors and dynamics, coordinated waters, and base-pair propeller twist angles[END_REF] and Class II Amynoacyl tRNA synthetase (Figure 2.7 right) [START_REF] Ruff | Class II aminoacyl transfer RNA synthetases: crystal structure of yeast aspartyl-tRNA synthetase complexed with tRNA (asp)[END_REF]. The sequences were taken from PDB database [START_REF] Westbrook | The Protein Data Bank and structural genomics[END_REF][START_REF] Berman | Announcing the worldwide protein data bank[END_REF]. Here, the section of the landscape is defined by two referent structures, which in this case is MFE and some high-energy structure that is distant from the MFE. The energy of secondary structures is indicated either by a color code, such as in the provided example, or by z coordinate in case of tri-dimensional visualization. The tri-dimensional variant then depicts all structures with higher energy than neighboring states as a peak, while locally optimal structures would be represented by a hole. The lowest point of the entire landscape would be MFE. Alternatively, the RNA folding landscape can be represented as a graph, with nodes being the structures and the edges representing the transitions. By definition, each transition is reversible.

The implication of the RNA folding landscapes in the RNA kinetics stems from the fact that besides the secondary structures themselves, it also represents their evolution, resp. folding pathways. A folding is a process by which a secondary structure forms or changes from one to another. Here, we see the folding as a hierarchical process where an RNA sequence changes its structure to one of the neighboring states. The folding pathway is then a path in RNA folding landscape consisting of secondary structures which RNA adopts when folding from one secondary structure to another. The objective of RNA kinetics is to study these The examples of a folding landscapes. The folding landscape for Yeast Phenylalanine tRNA (PDB ID 1TRA, left) [START_REF] Westhof | Restrained refinement of the monoclinic form of yeast phenylalanine transfer RNA. Temperature factors and dynamics, coordinated waters, and base-pair propeller twist angles[END_REF] and Class II Amynoacyl tRNA synthetase (1ASY, right) [START_REF] Ruff | Class II aminoacyl transfer RNA synthetases: crystal structure of yeast aspartyl-tRNA synthetase complexed with tRNA (asp)[END_REF]. The first and second reference are respectively MFE structure and a random high-energy structure distant from MFE structure. The picture was generated by an Rscript script that will be possibly available with the future releases of VIEN-NARNA library [START_REF] Lorenz | Vi-ennaRNA Package 2.0[END_REF].

folding pathways and the evolution of RNA structure, something for what the folding landscape is a great tool.

The methods to generate and study the RNA folding landscape can be divided into two categories:

• Step-by-step methods: The folding process is modeled as a Markov chain and simulated at a certain resolution, usually on base-pair [START_REF] Flamm | RNA folding at elementary step resolution[END_REF] or helix [START_REF] Ludmila V Danilova | RNAKinetics: a web server that models secondary structure kinetics of an elongating RNA[END_REF] level. However, since the number of the secondary structures raises exponentially with the length of sequence [START_REF] Zuker | RNA secondary structures and their prediction[END_REF] and the number of folding pathways also increases exponentially with it [START_REF] Flamm | Beyond energy minimization: approaches to the kinetic folding of RNA[END_REF], this type of analysis becomes quickly infeasible.

• Coarse-grained methods: The main idea is to simplify the folding landscape first. This necessitates to identify the key elementskey landmarks -of the folding landscape, the most important secondary structures. These are connected together in a manner to accurately represent the kinetics of a given RNA sequence. Such representation is then used to perform an analysis of RNA kinetics. Due to reduced number of structures and folding pathways, this strategy is applicable even to longer RNA sequences. However, its quality greatly depends on that of the first step, since the omission of certain key structures might have considerable impact on the analysis of kinetics and probably result in biased observations.

There is a number of currently existing approaches that generate these key landmarks [START_REF] Wuchty | Complete suboptimal folding of RNA and the stability of secondary structures[END_REF][START_REF] Li | Finding stable local optimal RNA secondary structures[END_REF][START_REF] Kucharik | Basin Hopping Graph: a computational framework to characterize RNA folding landscapes[END_REF][START_REF] William | Computing the partition function for kinetically trapped RNA secondary structures[END_REF] and assemble them into a folding landscape [START_REF] Flamm | Barrier trees of degenerate landscapes[END_REF][START_REF] Kucharik | Basin Hopping Graph: a computational framework to characterize RNA folding landscapes[END_REF]. Our work mainly concentrates on establishing a method capable of selecting the key landmarks of RNA folding landscapes in a more efficient manner than currently existing methods.

Selecting suboptimal secondary structures

For kinetic analysis to be possible for long RNAs, the folding landscape must be reduced to key suboptimal structures. The question is how to define and select them. This step is instrumental due to having the instrumental impact on the quality of the resulting landscape and adjacent analysis.

There are multiple strategies that we can chose from. The first one is to select only structures within a certain energy range from the MFE. Due to imprecisions of computations by Turner Energy Model, the computed MFE might not be MFE in reality, but it will not be too much far off. Second reason, the energy structures close to MFE are more stable and consequently a better candidates for a functional structure. However, there have been an exceptions to this, such as an artificially designed riboswitch where the energy difference between MFE and its metastable, functionally important, structure is almost 30 kcal.mol -1 [START_REF] Paul | RNA secondary structure: physical and computational aspects[END_REF]. Another thing to consider is that the number of structures selected this way still increases exponentially with n.

Another strategy is to only consider suboptimal RNA secondary structures S that are locally stable. The locally stable structures are considered more likely to be a candidates for an important functional structure due to the energy barriers that need to be surpassed for S to be changed, granting it a certain stability and allowing to sequence to exert its function. On the other hand even if the locally optimal secondary structure is not functionally important, it still can be crucial transitive state with the impact on the folding path and by extension on the result. For this reason, in the first part of this thesis we will work with locally optimal structures.

Unfortunately, like for the selection by energy band, the number of locally opti-mal structures also increases exponentially with n and since the folding landscape presents a high number of them, it can be problematic. Therefore, it is necessary to be able to select the best candidates among them. To achieve this, we resort to the sampling, which has for objective to select the most suitable locally optimal secondary structures according to a given criterium.

Boltzmann Distribution

The selection happens at two levels. The first is the selection of a subspace of S n (all secondary structures, only locally optimal ones etc.), while the second is to choose its best candidates. The second selection can be done either in deterministic or stochastic, random manner. Here we will concentrate on the stochastic sampling methods.

To perform a stochastic sampling, it is necessary to define a distribution over the entire space of secondary structures. The distribution is computed according to the value of S(S) by which we want the results to be filtered, such as the lowest free energy or biggest number of base pairs.

The distribution of secondary structures that is frequently used in sampling of RNA secondary structures is the Boltzmann distribution, sometimes also called Gibbs distribution. This distribution was formulated by Ludwig Boltzmann [START_REF] Sharp | Translation of Ludwig Boltzmann's paper "On the relationship between the second fundamental theorem of the mechanical theory of heat and probability calculations regarding the conditions for thermal equilibrium" sitzungberichte der kaiserlichen akademie der wissenschaften[END_REF] and later studied by Josias W. Gibbs [START_REF] Gibbs | Elementary principles in statistical mechanics[END_REF]. The Boltzmann distribution is used mainly in thermodynamics and is a probability distribution of various systems over different possible states. Since we consider states each aggregating of microscopic conformations of the molecules with the same secondary structure, and not separate conformations themselves, the Boltzmann distribution is indeed applicable on free energies.

Definition 2.5.4 (Boltzmann factor): et the system be the sequences w and the states the secondary structures S ∈ S n , each with free energy E S . The Boltzmann factor Z S Z S Z S of S is

Z S = e -E S k B T (2.9)
where T is an absolute temperature and k B is Boltzmann constant or perfect gas constant.

In our work, we consider k B ≈ 1.987 × 10 -3 kcal.mol -1 .K -1 .

Definition 2.5.5 (Partition function): partition function Z Z Z is a description of the statistical properties of the system, here RNA in ensemble of states from S n , in a thermodynamic equilibrium. It is computed as

Z = S∈Sn e -E S k B T . (2.10)
From this expression and the equation 2.9, we can obtain the probability of observing S.

Definition 2.5.6 (Boltzmann probability):

The probability p(S) p(S) p(S) of observing w folded in S ∈ S n under the assumption of Boltzmann distribution is:

p(S) = Z S Z = e -E S k B T X ∈Sn e -E X k B T (2.11)
These probabilities can be used to sample S in a random manner. Since S with low energy have higher Z S , they are more likely to be chosen, while the high energy structures still have non-zero probability to be selected since ∀S, Z S > 0. While priority is given to low-energy candidates, structures having higher free energy are still theoretically taken into consideration, and can be relatively favored by higher temperatures.

The first algorithm that computed the partition function for ensemble of secondary structures was John S. McCaskill [START_REF] Mccaskill | The equilibrium partition function and base pair binding probabilities for RNA secondary structure[END_REF]. He used the VZu algorithm (a non-ambiguous variant of Zuker algorithm, Equation 2.8) with transforming energy terms to partition function.

Completeness, Unambiguity and Acyclicity

Here we briefly mention three necessary properties of DP schemes that are used in algorithms computing the partition function [START_REF] Ponty | A combinatorial framework for designing (pseudoknotted) RNA algorithms[END_REF]. Both will be more detailed later. These three conditions are:

• Completeness: The space S n or its defined subsection is completely covered by the DP scheme, ie. there is not structure S that cannot be decomposed by said DP scheme while S ∈ S n .

• Unambiguity: Each secondary structure can be decomposed by the decomposition scheme in a single manner.

• Acyclicity: A DP scheme cannot present a decomposition where the element would be decomposed to itself and some other elements.

The DP scheme of Zuker and VZu algorithm presented here (Equation 2.8 and 2.7) completely covers S n -each nucleotide can be paired or not (as long as θ is respected), each loop can contain zero, one or more helices and each multibranch loop has to contain at least two helices possibly preceded and followed by unpaired bases. Note that the completeness depends on the definition of S n -the Zuker and VZu algorithms are complete when S n is defined as in Section 2.1, but not if pseudoknots are allowed.

In regards to unambiguity, the DP scheme of Zuker algorithm does not fulfill it, and is main reason why the second version was established. Suppose we have a segment WMi, j containing two helices Vi, k -1 and Vk + 1, j separated by a single unpaired nucleotide k. There are two ways how this local secondary structure can be decomposed:

WM i,j → WM i,k + V k+1,j → V i,k-1 + V k+1,j WM i,j → V i,k-1 + WM k,j → V i,k-1 + V k+1,j
This is not problematic when searching for MFE, since the minimum structure will be returned anyway, but when computing Partition function its Boltzmann factor will appear twice, introducing bias. For the VZu algorithm we have:

WM i,j = M i,k + M 1 k+1,j → M 1 i,k + C k+1,j → M 1 i,k-1 + C K+1,j → C i,k-1 + C k+1,j
Note that WM i,j is not present in VZu algorithm, it is just to point out the equivalence between the two versions. The second version respects unambiguity and can be used to compute the partition function, which is the base of McCaskill algorithm.

As for acyclicity, the presence of a cycle in DP scheme implies an infinite number of the possible decompositions. This makes their adaptation to scoring algorithms impossible. One possibility of verifying an acyclicity is therefore to see whether each decomposition can be done in finite number of steps. All of scoring algorithms presented in this work are acyclic.

McCaskill Algorithm

The partition function shows an interesting property with the respect to the decomposition of loop and the energies deduced from the Equation 2.5. Since E S is the sum of the energies of each of its loops E l , for Boltzmann factor we have:

Z S = e L∈S -E l k B T = l∈S e -E l k B T = l∈S Z l (2.12)
Consequently, the Boltzmann factor of a single structure can be decomposed into the product of contributions of each of the loops l. This can be used to transform the VZu algorithm into McCaskill algorithm. However, since in this case we want to compute the Partition function, we also need to sum the contribution of all secondary structures. We define:

• Z F i,j : Partition function over all S(i, j) ∈ S(i, j) • Z C i,j : Partition function over {S(i, j)|S(i, j) ∈ S(i, j), (i, j) ∈ S(i, j)}) • Z M i,j : Partition function over all S(i, j) ∈ S(i, j) where S(i, j) is a part of multibranch loop containing at least one helix

• Z M1

i,j : Partition function over all S(i, j) ∈ S(i, j) where S(i, j) is a part of multibranch loop containing exactly one helix These states are deduced from F i,j , respectively C i,j , M i,j , and M 1 i,j of VZu algorithm. The algorithm of McCaskill is then:

Z F i,j = Z F i,j-1 + i≤k≤j-θ-1 Z F i,k-1 × Z C k,j Z C i,j = e -E H (i,j) k B T + i<k<l<j k≤l-θ-1 e -E ILG (i,j,k,l) k B T × Z C k,l + i<k<j e -a-b k B T Z M i+1,k-1 × Z M1 k,j-1 Z M i,j = i<k<j e -c×(k-1) k B T Z M1 k,j + i<k<j Z M i,k-1 × Z M1 k,j Z M1 i,j = e -c k B T × Z M1 i,j-1 + e -b k B T × Z C i,j (2.13)
with E H (i, j), E ILG , a, b, c having the unchanged significance. The sums are explained by the need of summing up all contributions of each S. Since ∀w[i, j] | j ≤ i + θ the initial conditions are E(i, j) = 0, we have

∀{i, j | j ≤ i + θ}, Z i,j = Z B i,j = Z M i,j = Z M1 i,j = 1.
The value Z 1,n equals to Z. In the original paper [START_REF] Mccaskill | The equilibrium partition function and base pair binding probabilities for RNA secondary structure[END_REF], a partition function Z is instead denoted as Q, which has the same significance. The complexity is the same as for the VZu Algorithm -O(n 4 ) in time, which becomes O(n 3 ) if the maximum size of the internal loop becomes limited, and O(n 2 ) in space.

The algorithms computing Z of S n or its part like that of McCaskill are the first step towards the sampling of secondary structures.

Sampling and Stochastic backtrack

The sampling consists of the random selection of the secondary structures S from the space S n according to the probability distribution, in our case the Boltzmann distribution. The sampling algorithm employing McCaskill algorithm was first introduced by the works of Y. Ding and C.E. Lawrence [START_REF] Ding | A statistical sampling algorithm for RNA secondary structure prediction[END_REF]. It consists of two steps:

• Pre-computation: The computation of Z. More importantly, the matrices Z F , Z C , Z M and Z M1 are filled during this step, therefore ∀i, j, 1 ≤ i < j ≤ n, the values Z F , Z C , Z M and Z M1 are calculated.

• Stochastic backtrack: Performed once to generate a single S. The principle is similar to the backtrack as presented in Section 2.3.3, however instead of searching for a combination of a loop that gives a value given a matrix, we have a possibility to choose a certain loop l. Suppose we have w[i, j] and a possibility to choose a base pair, an unpaired base, a splitting point or a loop q 1 with a probability pair p(q 1 ) that occupies some bases and leaves sub-

strings w[i 1 , j 1 ], w[i 2 , j 2 ] • • • w[i x , j x ]
untouched. The value of p(q 1 ) is computed from DP matrices. If q 1 is chosen, then we can choose q 2 in the same manner for w[i 1 , j 1 ] with the probability p(q 2 ) and so on. The process is repeated until the substrings too short to contain a base pair. When no substring is left, the secondary structure is complete. The recursive application of these steps leads to a generation of S, from the search space S n , with Boltzmann probability.

The algorithm must give a way to access the probability values and these have to be stored in DP matrices. In the case of McCaskill algorithm, they can be obtained from Equation 2.13 where left-hand term represents the weight of all options for w[i, j] and the right-hand terms a list of choices available to be taken, giving access to p(q).

Since the introduction of sampling by Ding and Lawrence, other approaches have appeared, usually employing Turner Energy Model. These include a decomposition of a structure in stacks [START_REF] Li | Finding stable local optimal RNA secondary structures[END_REF] or by varying a temperature during sampling [START_REF] Kucharik | Basin Hopping Graph: a computational framework to characterize RNA folding landscapes[END_REF].

Formalization of Dynamic Programming Schemes

Since we just presented the state-of-the-art DP schemes employed for secondary structure prediction, it is a good moment to propose the formalism that would help us to generalize and unify them.

Definitions

Objective 2.6.1. Introduce a formalism that allows to describe DP schemes in an unified manner and allowing to explicitly manipulate the different computation steps. This formalism

• takes a DP scheme with all possible states and an entry state q root as an entry

• returns, for each state the associated accessible search space and scoring function expressed in an unified language that can be commonly used for DP schemes evaluation RNA secondary structures.

To do this, it is first necessary to identify the unitary actions of such DP scheme and define them in a broadest way possible. To avoid a confusion, we will, throughout this section, demonstrate how to apply it to formalize the Nussinov algorithm. The formalism presented here is inspired by and extends that introduced by Yann Ponty and Cédric Saule in relation to the representation of secondary structures as hypergraphs and pseudoknot analysis [START_REF] Ponty | A combinatorial framework for designing (pseudoknotted) RNA algorithms[END_REF].

A DP scheme relates the (optimal) score for a given subproblem, on the left-hand side of the equation, to the combination of scores achieved over the selected subset of subproblems on its right-hand side.

Definition 2.6.1 (DP scheme): A DP scheme is a tuple (Q, q root , ρ) such that:

• Q is the set of states q, eg. the couple of base pairs (i, j) associated with the right-hand-side terms of the scheme;

• q root ∈ Q is the initial state of the scheme, such as (1, n);

• ρ : Q → (Q * × D × R)
* is a function associating to each state q a set of derivations relating q to other states, with each derivation having attributed a constructor denoted λ, eg. an unpaired base, and a score to it.

Therefore a DP scheme relates the different states by a set of derivations. The notion of the derivation is also connected to that of constructors. We consecutively introduce all of these terms.

Definition 2.6.2 (State):

A given subproblem is, in proposed formalism, characterized by a state q q q ∈ Q, where Q denotes the set of all states with associated subproblems.

Therefore a given problem is modeled by a state q and its solution is indicated by a score over this state. The computation of a DP scheme total score requires a memorization of |Q| scores, pointing to an overall memory complexity of O(|Q|).

Example 2.6.1. For Nussinov algorithm, the states are any (i, j) such that 1 ≤ i < j ≤ n, where (i, j) is the couple of, not necessarily paired, nucleotides and we

have Q = {(i, j) | 1 ≤ i < j ≤ n}.
Each subproblem is completely defined by the nucleotide couple (i, j), therefore this definition of Nussinov states is sufficient.

Besides the score, each state q has a set of secondary structures s, that are accessible from this state, associated to it.

Definition 2.6.3 (Search space):

A search space is set of secondary structures s that are associated to a state q ∈ Q and can be accessed by the DP scheme from the sets of q 1 , . . . , q u ∈ Q. We denote it by L(q).

Note that L(q) does not contain only the structures that are merely visited, but all those that can be accessed from q. Also these structures do not have to have the length of n. For this reason we denote s ∈ L(q) a secondary structure belonging to a search space of q.

Example 2.6.2. A search space for Nussinov algorithm and sequence w, L(q) where q = (1, 5) contains all secondary structures that can be generated on w [START_REF] Akutsu | Dynamic programming algorithms for RNA secondary structure prediction with pseudoknots[END_REF][START_REF] Baker | Molecular biology of the gene[END_REF].

A structure S ∈ L(q) can be generated from some structures (s 1 , . . . , s u ) ∈ (L(q 1 ), . . . , L(q u ))

if q and q 1 , . . . , q u are related by the derivation.

Definition 2.6.4 (Constructor):

A constructor λ is a function λ : S k n → S n , k ∈ N, that generates a structure s s s from a vector of k (sub)structures s 1 , . . . , s k s 1 , . . . , s k s 1 , . . . , s k .

By abuse of notation, we might note by λ(X , Y) the construction λ(X, Y), with X ∈ X , Y ∈ Y. A constructor space is denoted D, the notation used in Definition 2.6.1.

Example 2.6.3. The Nussinov algorithm employs two different constructors: The constructor λ unp i,j where the secondary structure s is such that w[j] is unpaired and is applied to s ∈ L((i, j -1)), and the constructor λ pair i,j,k , where s is generated from substructures for states (i, k -1) and (k + 1, j -1) and (k, j) ∈ φ pair (w). In both cases i, j indicate the interval for which the constructor is employed, while in the second case k denotes the pairing partner to w[j]. The operator space is

D Nus = {λ unp i,j , λ pair i,j,k | 1 ≤ i ≤ k < j ≤ n}
The constructor is attributed to a specific derivation.

Definition 2.6.5 (Derivation):

A derivation q λ -→ c λ q 1 , . . . , q u is a function that associates a tuple ({q 1 , . . . , q u }, λ, c) to a state q, representing a derivation from q to q 1 • • • q u by λ and having cost c λ .

Therefore there are three layers to each derivation q λ -→ c λ {q 1 , . . . , q u }:

• Derivation q λ -→ c λ
{q 1 , . . . , q u } itself, relating the state q to q 1 , . . . , q u ;

• Construction of accessible search spaces, done by λ on {L(q 1 ), . . . , L(q u )};

• Computing the score of q, from scores of q 1 , . . . , q u and using the function employing c λ .

Note that some derivations might relate q to ∅, in which case we call them terminal derivations.

Each derivation q λ -→ c λ {q 1 , . . . , q u } corresponds to a specific case. In fact, there are two layers that can be distinguished for a given DP scheme:

• Cases: One DP scheme can include multiple cases on its right-hand side.

These cases correspond to the derivations. In the case of Nussinov algorithm, there is one case for w[j] being unpaired and a case for each pair (k, j) that can exist.

• Subproblems: Single derivation relates a single state q on the left-hand side to states q 1 , . . . , q u or subproblems. For λ unp i,j in the case of Nussinov algorithm, it connects one state on the left-hand side to one state on the right-hand side of the DP equation. In the case of λ pair i,j,k , it is one and two states respectively. Therefore, ρ(q) ρ(q) ρ(q) from the Defintion 2.6.1 is the set of derivations of form q λ -→ c λ {q 1 , . . . , q u } that relate q to q 1 , . . . , q u and generate s ∈ L(q) by constructor λ. The score of q is computed from the scores of q 1 , . . . , q u using value c λ . Alternatively we can also denote that ({q 1 , . . . , q u }, λ, c λ ) ∈ ρ(q). Example 2.6.4. We demonstrate that the DP scheme can be completely captured by the proposed formalism (Q, q root , ρ). Assume that Q, resp. D contain all q, resp. λ that appear in given DP scheme. Each DP equation relates a single state q ∈ Q on left-hand side and a number of different cases each relating a number of states q 1 , . . . , q u ∈ Q u on the right-hand side. Such case is completely defined by the derivation:

q λ -→ c λ {q 1 , . . . , q u }
where c λ is the energy bonus/penalty for using the constructor λ. The DP scheme is expressed by a union of all cases:

ρ(q) = q λ -→ c λ {q 1 ,...,q u } ({q 1 , . . . , q u }, λ, c λ )
Therefore as long as necessary λ and q are available, it is possible to define any derivation q λ -→ c λ {q 1 , . . . , q u } and all of its cases by ρ(q).

Example 2.6.5. In the case of the Nussinov algorithm, ρ(q) two following derivations are possible:

• (i, j) λ unp i,j ---→ 0 (i, j -1) -w[j] is unpaired; • q λ pair i,j,k ---→ 1 {(i, k -1), (k + 1, j -1)} where i ≤ k < j -w[j] is paired to w[k]
Consequently, the ρ Nus is:

ρ Nus ((i, j)) =                      ({(i, j -1)}, λ unp i,j , 0) if w[j] is unpaired i≤k<j (k,j)∈φ pair (w) {({(i, k -1), (k + 1, j -1)}, λ pair i,j,k , 1)}
Note that ρ(q) only contains the associations that are possible;

q λ pair i,j,k ----→ c λ pair i,j,k {(i, k -1), (k + 1, j -1)} does not exist if w[k]
and w[j] cannot be paired.

All possible structures s ∈ L(q) are generated by applying the constructor λ from derivation q λ -→ c λ {q 1 , . . . , q u } on (s 1 , . . . , s u ) ∈ (L(q 1 ), . . . , L(q u )).

Definition 2.6.6 (Generation of search space):

The search space L(q) can be generated by performing every possible derivation q λ -→ c λ {q 1 , . . . , q u }:

L(q) = q λ -→ c λ {q 1 ,...,q u } λ(L(q 1 ), . . . , L(q u ))
At any moment ∀q, L(q) ⊆ S.

With the search space L(q) and its construction defined, we have all elements to formalize the notion of completeness and unambiguity.

Definition 2.6.7 (Completeness): Let q root be an initial state and S n the space of all secondary structures of size n. The completeness implies:

L(q root ) = S n (2.14)
The completeness implies that every S ∈ S n is generated at one point or another. The space S n also might be restricted by specific conditions. In that case the completeness implies that the entirety of the restricted space is generated by DP scheme.

The unambiguity implies that every secondary structure is decomposed in an unique manner.

Definition 2.6.8 (Unambiguity):

Let x = q λ -→ c λ {q 1 , . . . , q u } and the L(x) = λ(L(q 1 ), . . . , L(q u )). A DP scheme is unambiguous that if x = q λ -→ c λ {q 1 , . . . , q u } and y = q λ --→ c λ
{q 1 , . . . , q v } are two derivations of the state

q such that x = y then L(x) ∩ L(y) = ∅.
In other words, said structure cannot be decomposed in multiple different ways, or the same s cannot be generated by two different λ 1 and λ 2 .

There is also an important condition of acyclicity of DP scheme. If DP scheme presents a cycle, it cannot be adapted for a scoring algorithm since it presents an infinite number of decompositions.

Definition 2.6.9 (Acyclicity): The DP scheme is acyclic if for any given sequence of derivations λ(L(q 1 ), . . . , L(q u )), each q is derived at most once.

In other words, no state can appear twice for the same decomposition.

Applications

Assume c(s) is the score of s. From ρ(q), we can define the function that computes the final score for every state q ∈ Q. As previously seen, there are two scoring strategies that are applied to DP schemes.

Problem 2.6.1:

• Input: DP scheme (Q, q root , ρ);

• Output: MinScore(q) such that:

MinScore(q) = min s∈L(q) (c(s)) .
The problem for maximization scoring is similar:

MaxScore(q) = max s∈L(q) (c(s)) .
Such scoring strategy is applied by Nussinov (maximization) [START_REF] Nussinov | Algorithms for loop matchings[END_REF] or Zuker (minimization) [START_REF] Zuker | Optimal computer folding of large RNA sequences using thermodynamics and auxiliary information[END_REF] algorithms. Here we present its global formalization.

With the scoring established, we can proceed to its formulate its computation using DP schemes.

Theorem 2.6.1 (Minimization/maximization scoring): Let (Q, q root , ρ) be a complete and unambiguous DP scheme. MinScore(q) can be computed in O(|Q| + |D|) time and O(|Q|) memory by the following recurrence:

MinS(q) = min q λ -→ c λ {q 1 ,...,q u } c λ + u i=1 MinS(q i ) (2.15)
The maximum score MaxS(q) is computed in O(|Q|+|D|) time and O(|Q|) memory by the following recurrence:

MaxS(q) = max q λ -→ c λ {q 1 ,...,q u } c λ + u i=1 MaxS(q i ) (2.16)
Proof 2.6.1. We demonstrate this expression by induction on the maximum number of derivations M that can be done on q.

Base case (M = 1): Consider a state q where all derivations take the form

q λ term ----→ c λ term ∅.
L(q) contains a set of structures s where each of them is generated uniquely by λ term on q. Consequently, the score c(s) is

c(s) = c λ term .
We assume completeness, L(q) is completely generated by λ term from all associated derivations described above. Therefore, the minimum for q is given by

MinS(q) = min q λ term ----→ c λ term ∅ c(s) = min s∈L(q) c(s).
It follows that MinS(q) = MinScore(q).

Induction (M = m): Assume that ∀i ∈ Z + , the state q i , can be derived at most m -1 times and that we have MinS(q i ) = MinScore(q i ).

Consider q and the derivation q λ -→ c λ {q 1 , . . . , q u }. Consequently, q can be derived at most m times. We have:

MinS(q) = min q λ -→ c λ {q 1 ,...,q u } c λ + u i=1 min s ∈L(q i ) c(s ) .
The term u i=1 min s ∈L(q i ) (c(s )) can be also interpreted as a lowest sum of free energies of the combination (s 1 , . . . , s u ) ∈ (L(q u ), . . . , L(q u )), which leads us to

u i=1 min s ∈L(q i ) c(s ) = min (s 1 ,...,s u )∈ (L(q u ),...,L(q u )) u i=1 c(s i ) .
Assume now we derive q by some derivation q λ --→ c λ {q 1 , . . . , q u }. We get

c λ + u i=1 min s ∈L(q i ) c(s ) = c λ + min (s 1 ,...,s u )∈ (L(q u ),...,L(q u )) u i=1 c(s i ) .
This term returns the minimal value corresponding to some combination of (s 1 , . . . , s u ) ∈ (L(q u ), . . . , L(q u )) and the contribution of c λ of derivation q λ --→ c λ {q 1 , . . . , q u }.

Therefore

c λ + min (s 1 ,...,s u )∈ (L(q u ),...,L(q u )) u i=1 c(s i ) = c(s λ )
where c(s λ ) = λ (s 1 , . . . , s u ) ∈ λ (L(q u ), . . . , L(q u )) that has minimum score along all secondary structures generated by λ , and s λ ∈ L(q). This applies to

every q λ -→ c λ {q 1 , . . . , q u }. Since c(s λ ) is minimal for q λ --→ c λ {q 1 , . . . , q u }, we can follow by MinS(q) = min q λ -→ c λ {q 1 ,...,q u } (c(s λ ))
There is another condition in DP scheme (Q, q root , ρ) being complete. In that case the application of every constructor λ generates L(q). Consequently, we can write

MinS(q) = min L(q) (c(s)) .
If follows that MinS(q) = MinScore(q). The demonstration for MinS(q) is done in the same way.

Example 2.6.6. Nussinov algorithm is the base pair maximization DP scheme.

The maximizing function can be split in cases depending whether w[j] is paired.

The maximum number of base pairs for S(i, j) ∈ L((i, j)) can be computed as MaxS((i, j)) = max( max

q λ pair i,k,j ---→ 1 {(i,k-1),(k+1,j-1)} (MaxS((i, k -1))+ + MaxS((k + 1, j -1) + 1), MaxS ((i, j -1)))
(2.17)

We can define the similar problem for the computation of the partition function.

Problem 2.6.2:

• Input: DP scheme (Q, q root , ρ);

• Output: PfScore(q) such that:

PfScore(q) = s∈L(q) e -c(s) k B T
.

As an example, the partition function is computed by McCaskill algorithm [START_REF] Mccaskill | The equilibrium partition function and base pair binding probabilities for RNA secondary structure[END_REF].

The formalization of partition function scoring of DP scheme (Q, q root , ρ) is similar as in the case of the score minimization/maximization.

Theorem 2.6.2 (Partition function scoring): Let (Q, q root , ρ) be a complete and unambiguous DP scheme. PfScore(q) is computed in O(|Q| + |D|) time and O(|Q|) memory by the following recurrence:

PfS(q) = q λ -→ c λ {q 1 ,...,q u } e -c λ k B T × u i=1 PfS(q i ) (2.18)
Proof 2.6.2. The demonstration can be made by induction on the maximum number M of successive derivations that can be performed on q.

Base case (M = 1): Consider any state q that can be derived maximum once. Therefore every derivation takes a form

q λ term ----→ c λ term ∅.
L(q) contains structures composed of an element generated by λ term . The Boltzmann factor of a single s is then:

Z s = e -c λ term k B T = e s k B T .
If we assume completeness, the set of s generated by all λ term gives L(q). Therefore

PfS(q) = q λ term ----→ c λ term ∅ e -c λ term k B T = s∈L(q) e -c(s) k B T .
since we assume unambiguity, meaning is no more derivation per the same structure than one. It follows that PfS(q) = PfScore(q).

Induction (M = m): Consider, ∀i ∈ Z + , the states q that can be derived maximum m -1 times and that they satisfy PfS(q ) = PfScore(q ). Assume now a state q derivation q λ -→ c λ {q 1 , . . . , q u }. This implies that q can be derived at most m times. We have

PfS(q) = q λ -→ c λ {q 1 ,...,q u }   e -c λ k B T × u i=1   s ∈L(q i ) e -c(s ) k B T     .
The term

u i=1   s ∈L(q i ) e -c(s ) k B T  
is also the cartesian product. This product contains the sum of all possible combinations of (s 1 , . . . , s u ) ∈ (L(q u ), . . . , L(q u )). Consequently

u i=1   s ∈L(q i ) e -c(s ) k B T   = (s 1 ,...,s u )∈ (L(q 1 ),...,L(q u )) u i=1 e -c(s i ) k B T .
Therefore, for a specific λ such that q λ --→ c λ {q 1 , . . . , q u }, we get:

e -c λ k B T × u i=1   s ∈L(q i ) e -c(s ) k B T   = e -c λ k B T × (s 1 ,...,s u )∈ (L(q 1 ),...,L(q u )) u i=1 e -c(s i ) k B T e -c λ k B T × u i=1   s ∈L(q i ) e -c(s ) k B T   = (s 1 ,...,s u )∈ (L(q 1 ),...,L(q u )) e -c λ - u i=0 c(s i ) k B T
The energy of s λ , s λ ∈ L(q) that was constructed by λ is equal to the free energy of each of its elements, here u i=0 c(s i ) and of the new element constructed by λ ,

c λ . Therefore -c λ - u i=0 c(s i ) = c(s λ ).
Consequently, we can write, for λ :

e -c λ k B T × u i=1   s ∈L(q i ) e -c(s ) k B T   = s λ =λ (s 1 ,...,s u )∈ λ (L(q 1 ),...,L(q u )) e -c(s λ ) k B T
This property holds true independently of q λ --→ c λ {q 1 , . . . , q u }. It remains to sum up the contribution for every q λ -→ c λ {q 1 , . . . , q u }:

PfS(q) = q λ -→ c λ {q 1 ,...,q u }      s λ =λ(s 1 ,...,s u )∈ λ(L(q 1 ),...,L(q u )) e -c(s λ ) k B T      .
If we assume the DP scheme (Q, q root , ρ) is complete and unambiguous, we can say that:

• ∀s ∈ L(q), the structure must be generated at least by one λ (completeness)

• ∀s ∈ L(q), the structure cannot be generated by more than one λ (unambiguity)

This implies that the set off all elements s λ ∈ λ(L(q 1 ), . . . , L(q u )) derived all possible derivations q λ -→ c λ {q 1 , . . . , q u } is equal to L(q), since each term of s appears exactly once in the given sum. It follows

PfS(q) = s∈L(q) e -c(s) k B T . (2.19)
If follows that PfS(q) = PfScore(q).

Formalizing VZu Algorithm

Here we work with the unambiguous version of the algorithm employed in the VIENNARNA package. In this case, defining the states solely by the couple (i, j) is insufficient, since the state also depends on whether (i, j) ∈ S(i, j) we are constructing or whether we are currently processing a multibranch loop branches.

For each (i, j), four different states q i,j indexed by the Equations 2.8 must be defined:

q F i,j = {(i, j)} q C i,j = {(i, j) | (i, j) ∈ φ pair (w)} q M i,j = {(i, j) | w[i, j
] is in multibranch loop containing at least one helix} q M1 i,j = {(i, j) | w[i, j] is in multibranch loop containing exactly one helix}

The set of all states in VZu Algorithm is then given as

Q Zuk = {q F i,j , q C i,j , q M i,j , q M1 i,j | 1 ≤ i < j ≤ n}.
To define D Zuk we need to create a specific derivation for each relation and case from the Equation 2.8. We employ both λ unp i,j and λ pair i,j,k . We need also to introduce a constructor for each of specific loops. Here we list all types of constructors for D VZu , considering S(i, j) ∈ S(i, j) is generated by it:

• λ unp i,j : w[j] is unpaired; • λ pair i,j,k : (k, j) ∈ S(i, j); • λ H i,j : w[k] is unpaired ∀i < k < j, S(i, j) = {(i, j)}; • λ ILG i,j,k,l : , {(i, j), (k, l)} ⊆ S(i, j)) with i < k < l < j; • λ ML i,j,k : (i, j) ∈ S(i, j), at least one helix on w[i + 1, k -1], exactly one helix on w[k, j -1]. • λ ML2hel i,j,k : w[i, j] in multibranch loop, at least one helix on w[i, k -1], exactly one helix on w[k, j]. • λ ML1hel i,j,k : w[i, j] in multibranch loop, {(i , j )} ⊆ S(i, j) with i ≤< k ≤ i < j ≤ j, {w[l]i ≤ l ≤ k} is unpaired; • λ p i,j : w[i, j] in multibranch loop (i, j)S(i, j)
The D VZu is then defined as an union of these constructors for all possible combinations of corresponding parameters.

Only some of the constructors can be used to generate L(q). The relations are defined by DP scheme used in VZu algorithm. From it we can define ρ:

ρ VZu (q F i,j ) =          (q F i,j-1 , λ unp i,j , 0) i≤k<j (k,j)∈S(i,j) ({q F i,k-1 , q C k+1,j-1 }, λ pair i,j,k , 0) ρ VZu (q C i,j ) =                    (∅, λ H i,j , E H (i, j)) i<k<l<j (k,l)∈S(i,j) (q C k,l , λ ILG i,j,k,l , E ILG (i, j, k, l)) i<k<j ({q M i+1,k-1 , q M1 k,j-1 }, λ ML i,k,j , a + b) ρ VZu (q M i,j ) =          i<k<j ({q M1 k,j }, λ ML1hel i,j,k , c × (k -i)) i<k<j ({q M i,k-1 , q M1 k,j }, λ ML2hel i,j,k , 0) ρ VZu (q M1 i,j ) =    ({q M1 i,j-1 }, λ unp i,j , c) ({q C i,j }, λ p i,j , b) (2.20)
With ρ VZu defined, we complete the formalization of DP scheme used in VZu algorithm by the triplet (Q Zuk , q F i,j , ρ Zuk ). Here, q root = q F 1,n , and if completeness is validated, it stores the free energy of MFE structure. The recursion function to compute MFE can be obtained from this scheme by applying the Theorem 2.6.1 to it. Such a formalization can be used to establish the overall complexity of the algorithm. The complexity of an energy minimization algorithm is O(|Q| + |D|) in time and O(|Q|) in space (see Theorem 2.6.2). In the case of time complexity, the determining factor is the number of constructors for an internal loop, due to con-taining four parameters with each taking approximately n values, the most of all constructors, and more than |Q|, which is

O(n 2 ). Therefore O(|Q VZu | + |D V Zu|) = O(n 4
). This can be reduced to O(n 3 ) if the maximum number of unpaired bases in an internal loop/bulge is limited, because in that case one of the parameters can be expressed in function of the limit and the other three parameters, and for other constructors their number is majored by n 3 . The memory complexity is O(n 2 ), since all states can be stored in a limited number of two-dimensional matrices, therefore O(n 2 ).

Formalizing McCaskill Algorithm

The McCaskill Algorithm is similar to the VZu algorithm, the main difference is that the minimizing function MinS is replaced by the function computing the partition function PfS. Besides that, Q MC = Q VZu and D MC = D VZu and ρ MC = DVZu, therefore VZu and McCaskill algorithm are formalized in the exactly same way. These formalizations can be easily employed to prove the unambiguity of both algorithms, though it is required only in the latter case.

Example 2.6.7. We will demonstrate the unambiguity on VZu algorithm. The unambiguity implies a single s ∈ S is not generated by a two different constructors. This must be proven separately for each state q. We simply look what constraints different λ impose on implicated nucleotides.

In the case of L(q F i,j ), the proof is straightforward since λ unp i,j (L(q F i,j-1 )) implies that w[j] is unpaired while λ pair i,j,k (L(q F i,k-1 ), L(q F k+1,j-1 )) necessitates the presence of pair (k, j), and since base triplet is forbidden by definition, we have, for any i, j, k 1 , k 2 such that 1 < i ≤ k 1 < k 2 < j ≤ n:

λ unp i,j (L(q F i,j-1 )) ∪ λ pair i,j,k 1 (L(q F i,k 1 -1 ), L(q F k 1 +1,j-1 ) = ∅ λ pair i,j,k 1 (L(q F i,k 1 -1 ), L(q F k 1 +1,j-1 ) ∪ λ pair i,j,k 2 (L(q F i,k 2 -1 ), L(q F k 2 +1,j-1 ) = ∅
In a similar manner, each of the cases computing L(q C i,j ) is mutually exclusive; λ HP i,j , λ ILG i,j,k,l , λ ML i,j,k imply respectively zero, one or multiple base helices for w[i + 1, j -1] and each time the values of eventual coefficients k and l are different.

Note that in the last case the construction of q M1 i,j implies that the last helix in given multibranch loop begins at w[k], implying the mutual exclusivity between different structures generated by λ ML i,j,k for different values of k.

In the case of generating L(q M i,j ) the segment w[i, k -1] is either completely unpaired, or hosts another base pair. Like in previous case w[k] indicates the beginning of last helix in multibranch loop, therefore the cases cannot generate the same structure for different k, and λ ML1hel i,j,k and λ ML2hel i,j,k cannot result in same S(i, j) either due to difference in w[i, k -1].

The explanation for L(q M1 i,j ) is similar to the one written for L(q F i,j ), the difference being only base pair (i, j) being considered. This implies that for no state q two different constructors λ can generate the same S, therefore VZu algorithm is unambiguous. The same applies to McCaskill algorithm due to it using the same decomposition.

Stochastic backtrack

The general idea behind the stochastic backtrack was explained in Section 2.5.6. With the introduced formalism the central principle can be explained in a much more simplified and compact manner. The stochastic backtrack can be only performed by a DP scheme (Q, q root , ρ) used to compute the partition function, ie. employs the scoring function PfS. For each state q, we choose what derivation will take place.

Theorem 2.6.3 (Backtrack probability): By using the probability p(λ | q) of performing a derivation q λ -→ c λ {q 1 , . . . , q u } when currently being in state q and whose value is

p(λ | q) = e -c λ k B T × u i=1 PfS(q i ) PfS(q) (2.21)
the stochastic proces ulitmately terminates and returns an object s = λ(Lq 1 , . . . , Lq u ) with boltzman probability within L(q). Proof 2.6.3. From Proof 2.6.2 we know that contribution of a single derivation q

λ -→ c λ {q 1 , . . . , q u } towards PfS(q) is e -c λ k B T × u i=1   s ∈L(q i ) e -c(s ) k B T   .
It follows directly from Theorem 2.6.5 that PfS(q) gives the sum of contributions of all λ and that the probability p(λ | q) is computed as shown in Equation 2.21.

The stochastic backtrack of an entire structure is then performed by starting at q root and repeating the process for every state q 1 , q 2 , . . . , q u that gets selected.

Note that the order of treatment of states does not matter as long as it is deterministic, but they must be all treated.

The relation between RNA thermodynamics and RNA kinetics

We mentioned that due to the structural diversity of most larger RNA sequences, the RNA folding landscape is obtained by sampling from a representative set of structures S ∈ S n or some specific subspace as locally optimal structures, free energies computed using Turner Energy Model. This is followed by the study of folding pathways and finally by the analysis of the evolution of the entire system. Therefore, we use a thermodynamic approach to study and simulate the kinetics.

It is worth mentioning the relation between thermodynamics and kinetics. The folding of an RNA can be seen as a type of chemical reaction with the secondary structures instead of compounds. Here we will consider a following transformation of secondary structure S A of an RNA sequence into the state S B , and we assume its kinetics follows the same stochastic model as a kinetics of chemical reactions [START_REF] Flamm | RNA folding at elementary step resolution[END_REF]:

S A → S B
While thermodynamics, here represented by the parameters of Turner Energy Model, will allow us to know about the stability of the different secondary structures and consequently whether the reaction should be possible but nothing about its speed, the kinetics indicates how fast the transformation would happen, but says nothing about the equilibrium state. The thermodynamics therefore helps us to choose the possible transformations, which we can then investigate from the perspective of the kinetics.

However, since we compute the free energy using thermodynamic energy model, the obtained energies, and notably associated Boltzmann distribution are correct while in the state of thermodynamic equilibrium, but not outside of it. If the system still evolves, then the thermodynamic energy model does not apply.

Example 2.7.1. Suppose a system where all RNA molecules are unfolded at some starting time t start . According to Boltzmann distribution their proportion should be with few exception very small or zero. However, since we are obviously not in equilibrium, their proportion is 1.

On the other hand, the Boltzmann distribution can still provide a good approximation on the kinetics because it indicates towards what state the system evolves.

The main problem with the use of thermodynamic model for the study of kinetics is the lack of so-called activation energies E a . It can be, however, computed from the energies of the structures participating in transformation. Its values do not depend only on the energy of those structures between which the transformation occurs, but also whether they are neighbors or not. Consider the reaction between two neighboring structures

S A → S B ,
The activation energy of neighboring structures the E a is the energy that is necessary to supply in order to transform S A to S B and its value is:

E a = E S B -E S A (2.22)
Activation energy is sometimes also denoted as an energy barrier. The case of neighboring structures is simple to compute, because there is no intermediate step. If the activation energy is negative, then the transformation is spontaneous -happens on its own -and the energy is instead freed.

In the case where S A and S B are not neighbors, it depends whether the highest free energy of intermediate product S peak from transformation of S A to S B is higher than that of S B or not.

Definition 2.7.1 (Activation energy): Suppose a transformation between states

S A → S peak → S B
where S peak is the intermediate structure with the highest free energy and S peak > max(S A , S B ). The activation energy necessary for this transformation to overcome the energy gradient is

E a = E S peak -E S BA . (2.23)
If there are multiple folding pathways between distant S A and S B , the one considered is that with the lowest E S peak , since it is the most optimal. This way we compute an approximation of the E a .

Such approach has an advantage of being easily computable. The Turner Energy Model is quite complete and contains extensive tabulated values as well as relations to compute the free energy of different substructures. The free energy is relatively easy and fast to compute, an important factor when determining the energies of longer sequences. While the computation of the activation energy for the most optimal folding pathway between S A and S B is an NP-hard problem [START_REF] Ma Ňuch | NPcompleteness of the energy barrier problem without pseudoknots and temporary arcs[END_REF], it belongs to PTAS class (polynomial-time approximation schemes exist for it) [START_REF] Rafiey | Ordering with precedence constraints and budget minimization[END_REF] and a heuristic methods allowing to compute their approximate values exist, such as findpath from VIENNARNA or RNATABUPATH [START_REF] Dotu | Computing folding pathways between RNA secondary structures[END_REF]. The partition function and associated probabilities can be also computed in a polynomial time using DPbased programming. There are more precise methods such as molecular dynamic simulations [START_REF] David R Bell | Capturing RNA folding free energy with coarse-grained molecular dynamics simulations[END_REF], but they are difficult to implement and notably very demanding in terms of time computation, consequently being unusable for the analysis of longer RNA sequences. Also, even the approximate analysis using thermodynamic model can give us deeper insight on the RNA kinetics.

The thermodynamic model also describes the ideal state of the system, therefore the state towards which it constantly evolves and which it would achieve if given enough time. This alone can give us indication which folding pathways could be crucial, and complemented by the computation of approximate activating energies we can obtain sufficient estimation of the kinetics of the system. This, along with the relative simpleness of the computations, is the reason why we will use thermodynamic Turner Energy Model to study the kinetics of RNA folding.

Kinetics of RNA folding landscape

The RNA landscape provides a great tool to study RNA kinetics. Since the size of S n makes building and analyzing the RNA folding landscape impossible to perform in reasonable time, the first objective when studying an RNA folding landscape is to simplify it to only the key elements. This applies to the folding pathways as well. Here we explain the general approach to reduce RNA folding landscape which will be followed through this work.

Simplification of a landscape

One way to simplify the folding landscape is to represent only its locally optimal structures Therefore, here we consider the space of only locally optimal secondary structures. Let S n L ∈ S n be the space of all locally optimal secondary structures S L of length n. Each S L is connected to a number of secondary structures with higher free energies such that the free energy from said structure to S L is strictly descending along the path.

E S peak = min S∈ B S LA ∩ B S LB (E S ).
Such saddle point is also called direct saddle point. An indirect saddle point is a saddle point that can be obtained for two basins by passing through another basins. Note that the energy of indirect saddle of two basins may be lower that that of the direct saddle. This greatly impacts what path will be taken when transformation S n L A → S n L B takes place.

The activation energy associated to transformation

S n L A → S n L B is E a = E S peak -E Sn L A .
The energy of indirect saddle is obtained by taking the highest energy of S peak that is met along the path from basin to basin.

Since the saddles are not locally optimal structures, they are not included in S n L , yet they are crucial to the estimation of the kinetics. However, many methods determine the saddle structures in order to estimate E a . An example of such software is findpath from VIENNARNA or RNAtabupath [START_REF] Dotu | Computing folding pathways between RNA secondary structures[END_REF].

Since the folding landscape has a ragged nature, it presents many of small holes and therefore many of S L . This number would still be prohibitively large for the analysis. This is why the sampling methods are employed, with Boltzmann sampling studied in this work prioritizing lower free energy structures that are potentially more interesting.

Alternatively, it is possible to investigate the raggedness of the landscape beforehand by dedicated methods and use it to local minima selection. One possible approach is to compute persistence of each local minimum -the lowest activation energy that allows to access the minimum with even lower energy [START_REF] Cazals | Conformational ensembles and sampled energy landscapes: Analysis and comparison[END_REF][START_REF] Carr | Energy landscapes and persistent minima[END_REF]. Such values can be then used to select the local minima with high persistence to obtain their representative set. This method is mentioned for the sake of completeness and it is not further developed here.

The reduction of the number of structures used to model a folding landscape also diminishes the number of the folding pathways. Plus, since we consider only one the most optimal path between S LA and S LB where B S LA and B S LB are neighbors, the number of folding pathways is reduced even further.

However, it is not always known what basins are neighboring. For most of the existing software, this is not too much of problem since they compute the activation energy between two states without having to know that, but computing energy barrier is computationally rather demanding, therefore we cannot compute the energy barriers between every two existing states.

C(S LA ) = {S L ∈ S n L | d(S LA , S LB ) < c d , S L = S LA }
The E a are then computed for each C(S L ).

The c d as well as d function must be defined in a manner that the resulting landscape is connected, otherwise it will be impossible to perform a simulation due to the impossibility to access a part of the landscape from a specific structure.

The final result is the model of the RNA folding landscape reduced to locally optima structures and a key transitive structures with a corresponding set of folding pathways. While this model might be overly simplified, it allows already to study the effect of locally optimal structures on the kinetics of a given RNA sequence, something we are primarily interested in this work. The advantage is that the analysis of such landscape is faster than more complex models. For this reason this was the method we employed during the course of our study.

RNA Kinetics Study

Here we consider S ∈ S n L . The study may be performed on structures from S n in general, however in our context considering only S n L is sufficient. For each S L , we define P(S L , t), the proportion of RNA molecules of the system folded in S at the time t. We have ∀t,

S X ∈Sn L P(S X , t) = 1.
The folding pathways between the different S L can be seen as an iterative transformation implicating a certain number of intermediate structures [START_REF] Flamm | RNA folding at elementary step resolution[END_REF]. In the absence of inertia, the entire process can be seen as memoryless, since the choice of the next structure does not depend on the previous choices. Therefore, the changes of the proportions within the system can be modeled by a continuous time discrete state Markov Chain process. If we suppose the transitions between S LA and C(S LA ) and all transitions are reversible, then each P(S LA , t) depends on all influxes from C(S LA ) and all outfluxes towards this space. The evolution of P(S LA , t) can be captured by a chemical master equation (CME) [START_REF] Lecca | Stochastic chemical kinetics[END_REF].

dP(S L A , t) dt = S Ln ∈C(S LA ) P(S Ln , t) × k S Ln →S LA -P(S LA , t) × k S LA →S Ln . (2.24)
The values k S Ln →S LA are called transition rates. These describe the speed at which one structure is transformed to another. For the self-transition rate of S LA (rate for the S LA not being changed), we have [START_REF] Timm | Random transition-rate matrices for the master equation[END_REF]:

k S LA →S LA = - S Ln ∈C(S LA ) k S Ln →S LA (2.25)
If the system is in equilibrium at t → +∞, then dP(S LA , t) = 0 since system does not evolve, and also P(S L , t) = p(S), therefore P(S L , t) can be computed using Equation 2.11. Assuming equilibrium, we can utilize CME to compute the ratio of the transition rates for the same transformation in opposite directions:

e -E S Ln k B T × k S Ln →S LA Z = e -E S LA k B T × k S LA →Sn Z k S LA →S Ln k S Ln →S LA = e -E Sn k B T e -E S LA k B T k S LA →S Ln k S Ln →S LA = e -∆E S Ln -S LA k B T
with E S Ln , E S LA the energies of S Ln and S LA respectively and ∆E Sn-S A is their difference. Consequently, the ratio of the transition rates is the Boltzmann factor of their energy difference.

The computation of the transition rates is the necessary step before performing a simulation of a kinetics of specific RNA. The entire folding process is modeled by Markov Chains and the transition rates are used to establish a transition rate matrix that is necessary to perform a numerical integration of the Markov chains.

The simplest way to approximate their values is to use the Metropolis rule. In general, this rule assumes that that the transformation is guaranteed if the energy of the final state is lower than the energy of starting state:

k Sx→Sy =    τ × e -Ea x→y k B T if E ax→y > 0 τ otherwise (2.26)
where E ax→y is the activation energy of transformation S x → S y and k Sx→Sy the associated transition rate.

If S Lx , S Ly ∈ S 2 n L , the E ax→y ≥ 0, otherwise one of S L is not locally optimal. Consequently, for the transformation S Lx → S Ly of locally optimal structures the Metropolis rule can be reduced

k S Lx →S Ly = τ × e -Ea x→y k B T
In the case of locally optimal structures, this expression is in fact similar to Arrhenius law [START_REF] Arrhenius | Über die dissociationswärme und den einfluss der temperatur auf den dissociationsgrad der elektrolyte[END_REF], another law that can be used to compute the transition rates:

k S Lx →S Ly = A × e -Ea x→y k B T (2.27)
The last rule that can be and frequently is used to compute the transition rate is Kawasaki rule:

k S Lx →S Ly = τ × e -Ea x→y 2 boltzT
(2.28)

During this work, the energy barriers were computed using findpath for the pairs of locally optimal secondary structures. The details are specified in Section 3.3.4.

Coarse-grained model of a folding landscape

The different approach to simplify an RNA folding landscape and compute the transition rates is to introduce the macrostates. These macrostates represent multiple structures instead of the single one [START_REF] Michael | Efficient computation of RNA folding dynamics[END_REF].

We might consider each basin of attraction as a singular macrostate, reducing the number of states as well as that of folding pathways. In that case the macrostate associated to the basin B S L represents S ∈ B S L . The partition function of such macrostate is

Z B S L = S∈B S L Z S .
It follows the computing the free energy of the ensemble E B S L of the basin by:

E B S L = k B T ln(Z B S L ).
The transition rate between two macrostates k B S LA → k B S LB computes a bit differently than between two structures [START_REF] Michael | Efficient computation of RNA folding dynamics[END_REF]. If we know the transition rates k Sx→Sy such that S x ∈ B S LA and S y ∈ B S LB then

k B S LA → k B S LB = Sx∈B S LA Sy∈B S LB k Sx→Sy Z Sx Z B S LA
There is a number of different software that is capable of generating a model of an RNA landscape from the set of secondary structures along with the computation of the transition rate matrix using here described approach. They naturally implement a heuristic method to estimate the activation energies.

The most notable is barriers [START_REF] Flamm | Barrier trees of degenerate landscapes[END_REF], that creates a barrier tree. This tree is constructed from locally optimal structures by flooding algorithm -the locally optimal structures are the leaves of the barrier tree. Whenever two basins associated to the locally optimal structures merge a saddle point is found and a node connecting those two branches is placed. This is repeated until the root of the tree is found. The barrier tree graphs then represents the activation energies or energy barriers between two structures by the vertical length of the edges, and the free energy of locally optimal structures and saddles is determined by their vertical position.

Another approach to obtain a model of energy landscape is to compute the persistence diagram [START_REF] Cazals | Conformational ensembles and sampled energy landscapes: Analysis and comparison[END_REF][START_REF] Carr | Energy landscapes and persistent minima[END_REF]. This type of diagram depicts each local minimum and its connection to the local minimum with lower energy and minimum activation energy -the activation energy that effectively determines the persistence. Such visualization is quite compact and allows for fast identification of the most persistent local minima.

An analysis, more precise than that by barriers software, is performed by Basin Hopping Graph (BHG) [START_REF] Kucharik | Basin Hopping Graph: a computational framework to characterize RNA folding landscapes[END_REF]. This graph connects all neighboring basins B, represented by nodes and the edges represent the transitions along with the corresponding value of the activation energy E a . Unlike the barrier tree, which shows only the saddle points between two closest basins on their fusion, the BHG shows all transition paths and therefore does not succumb to a loss of an information. On the other hand its disadvantage is that it may return multiple additional transitive structures generated by heuristics detecting the saddle points, which will make kinetic simulation last longer.

This approach is more precise than just using locally optimal structures but is also much more demanding due to potentially introducing the number of new supplementary structures, which makes the analysis of multiple landscapes or landscapes of longer sequences time consuming.

Kinetic Simulation

The kinetic simulation can be performed once the transition rate matrix was established. This matrix has rows and columns represented by a certain secondary structure that is present within the landscape. The transition rates between the secondary structures that are not directly connected are supposed to be zero. The numerical integration is usually done by using the software dedicated for this purpose. During this work, we used treekin [START_REF] Michael | Efficient computation of RNA folding dynamics[END_REF].

The method employed in treekin is specifically the Metropolis-Hastings (MH) algorithm [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF], one of Markov Chain Monte Carlo methods. The stochastic algorithm mostly employed to determine a non-standard distribution, it can be also used, due to its iterative principle, to follow an evolution of the system towards its equilibrium as well. However, to employ the MH algorithm, the transition rates must fulfill the condition of detailed balance. Applying the usual, rejectionbased MH algorithm does not always fulfill this condition [START_REF] Flamm | Beyond energy minimization: approaches to the kinetic folding of RNA[END_REF]. However, the algorithm can modified into a rejection-less variant that is similar to Gillespie algorithm employed in chemistry [START_REF] Daniel | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF].

Assume there is a sequence w which can fold into a certain number of S L ∈ S n L .

The simplified general principle is as follows

• Initialization: For t 0 , have w folded in structure S LA .

• Iteration: Perform the following:

-Find all structures S Ln ∈ C(S LA ). Compute

k T = S Ln ∈C(S LA ) k S LA →S Ln .
-Select one from them, S LB , with the probability

p(S LB ) = k S LA →S Ln k T
-Make S LB the new structure of w.

-Move to a new time t 1 = t 0 + ∆t such that ∆t is randomly generated from an exponential distribution with the rate parameter k T .

-Repeat iteration steps until some specified t final is reached. This is then generalized on all secondary structures in the system in specific initial distribution. The simulation is then performed on the entire system.

The choice of the initial distribution of structures depends on the objective of the kinetic study. We may just want to study the general evolution of the system. For riboswitches, we may want to observe how the concentrations of the MFE and the metastable state. The initial distribution must be adjusted for these specific objectives.

Sampling Caveats

While the field of RNA kinetics studies rapidly advanced in few last years, there is still much to be done. Kinetic analysis of especially long sequences still represents a problem that is mainly related to the richness of the space of the secondary structures.

The Boltzmann sampling prioritizes the secondary structures with lower free energy. This probability decreases exponentially as the energy increases. Let -∆G be a free energy difference between S 2 and S 1 . We call by a frequency ratio η(S 1 , S 2 ) a ratio of probabilities p(S 1 ) and p(S 2 ) of observing secondary structures S 1 and S 2 respectively

η(S 1 , S 2 ) = e -E S 1 k B T e -E S 2 k B T = e E S 2 -E S 1 k B T = e -∆G k B T . (2.29)
While it is desired to sample low energy structures in priority, it may become problematic to generate a higher free energy structures if we want to do so. For ∆G = 1, η(S 1 , S 2 ) = 5.067, meaning per each 5 sampled S 1 , we observe one molecule in S 2 . For ∆ = 10 this value increases to 11147271. It quickly becomes impossible to access high energy structures without extensive sampling and even then it depends on the random number generator whether such structure will be selected.

A number of sampling strategies that tried to address this problem. As already stated, one method is to ample only locally optimal secondary structures S n L . This can be checked regularly at each step of decomposition. For example, the software RNAlocopt [START_REF] Lorenz | Vi-ennaRNA Package 2.0[END_REF] verifies whether the addition of a base pair will not result in the structure not being locally optimal anymore. RNASLOpt [START_REF] Li | Finding stable local optimal RNA secondary structures[END_REF] on the other hand assumes the stability is due to helical regions, which the authors call a stack, and searches for the stacks of certain length. For obtained S L , their values of E a are checked against a user-defined threshold whether they are stable enough.

However, while S n L is substantially smaller than S n , it still raises exponentially with n [START_REF] William | Computing the partition function for kinetically trapped RNA secondary structures[END_REF], therefore restricting S n to S n L will not completely solve the problem.

Another way to tackle this problem is to reduce the redundancy between of the sampled structures -sample the same structure repeatedly. One such approach is employed in software RNAlocmin [START_REF] Kucharik | Basin Hopping Graph: a computational framework to characterize RNA folding landscapes[END_REF], which uses the ξ-scheduling. The core of this approach is to sample multiple times on different temperatures, computed by ξ-scheduling, without scaling the Turner Energy Model-related factors -only the temperature T of Boltzmann factor is scaled. This results in smaller differences between the probabilities of different S and for temperatures sufficiently high temperature the Boltzmann distribution becomes uniform since the exponential factor asymptotically reaches 0. RNAlocmin exploits this to sample S at gradually higher temperatures, since smaller differences between probabilities mean higher probability to sample not-yet seen secondary structure. On the other hand the structures estimated at higher temperature are likely to have their energy overestimated and the redundancy is still present to some degree.

The way to prevent the repeats of samples altogether is to introduce the concept of non-redundant sampling. One of main objectives of this work is to research a sampling method that allows to generate each sample at most once without introducing bias in the distribution of not yet sampled structures. In next chapter, we will present the details on the concept of the method itself as well as how to integrate it in most efficient and easiest way in multiple algorithms. The efficient implementation is particularly important in regards to reduce the time consumed by sampling. We will mainly describe its implementation to Saffarian algorithm [START_REF] Saffarian | RNA locally optimal secondary structures[END_REF], a state-of-the-art algorithm that will be detailed in next chapter as well.

67

Chapter 3

RNANR: An Algorithm for Non-Redundant Sampling of RNA Secondary Structures

The main objective of this work is to establish an efficient method to study the kinetics of RNA sequences. We decided to achieve this by researching a nonredundant sampling method for the purpose of RNA secondary structure prediction. Since the quality of samples greatly impacts the quality of the RNA folding landscape, the main tool used to study RNA kinetics, it is reasonable to concentrate on amelioration of this important step. In plus, the non-redundant sampling method should be as optimal as possible, easy to implement, therefore easily separable form DP scheme it is adjacent to, and compatible with most DP schemes that compute partition function.

The algorithm that was used to research the non-redundant sampling is the algorithm based on combinatorics developed by Azadeh Saffarian, Hélène Touzet et al [START_REF] Saffarian | RNA locally optimal secondary structures[END_REF]. We will reference this algorithm as Saffarian algorithm. However, the algorithm itself does not compute the partition function nor includes the Turner Energy Model terms. After thoroughly introducing the mentioned algorithm and its formalization using the notions from Section 2.6, we will explain and validate its necessary modifications to compute these values and most importantly the probabilities p(S). We will then follow with explaining of the concept of nonredundant sampling.

The non-redundant sampling necessitates two things to function correctly. First, the distribution of not-yet sampled S after forbidding a sampled one must stay unbiased. This means the frequency ratio η(S 1 , S2) of any not already sampled S 1 and S 2 must stay same for each sampling iteration. Second, it must track all samples that were already generated. This also includes the chain of decision, or constructors λ, that were used to generate said structure. For this purpose, a specific data-structure that must be developed and implemented into the DP scheme to enable the non-redundant sampling. We will develop on both points in detail.

We mentioned that this approach should be compatible with most algorithms computing partition function and used for RNA structure prediction. We will support this point by demonstrating how one can implement the non-redundant sampling into RNAsubopt from VIENNARNA based on McCaskill algorithm.

Saffarian Algorithm

The Saffarian Algorithm, developed by Azadeh Saffarian, Hélène Touzet et al [START_REF] Saffarian | RNA locally optimal secondary structures[END_REF], is a combinatorial algorithm that returns only the locally optimal secondary structures. This algorithm does not compute the energy of any locally optimal secondary structure, instead it makes an assumption on local optimality based on combinatorial principle and then returns all the structures that satisfy that specific condition. The base version of the algorithm performs an exhaustive enumeration of all possible locally optimal secondary structures.

State of the Art and the general principle

The central hypothesis upon which the Saffarian algorithm is based is following: the structure is locally optimal only if it is saturated. In other words, there is no possible way to add a new base pair to the saturated structure S SAT . This however does not mean that S SAT contains the maximum base pairs as defined by Nussinov Algorithm, just that there is no base pair that can be added to the base pair configuration of S SAT .

The structures S n SAT delimit their own space.

Definition 3.1.2 (Saturated structure space): The saturated structure space is denoted by S n SAT ⊆ S n . The space of the saturated secondary structures S SAT (i, j) covering w[i, j] is denoted by S SAT (i, j). The space of all saturated secondary structures is denoted by S SAT with

S SAT = 1≤i<j≤n S SAT (i, j).
The reason that the saturated structures can be considered as locally optimal is that the base pairs have stabilizing effect on the structure. However, this is not always true, especially for the lonely base pairs -the base pairs that are not stacked with another base pairs. In fact, most of the stability of the structure comes from the helices, or stacked bas pairs [START_REF] Paul | RNA secondary structure: physical and computational aspects[END_REF]. To take this factor into the account, the Saffarian algorithm allows only for helices H(i, j, lh) with a certain minimum length a. In the rest of this chapter, we suppose a = 3.

If the minimum helix length is imposed, the definition of the saturated structure is slightly redefined: Definition 3.1.3 (Parametrized saturated structure): If all S can contain only helices H(i, j, lh), 1 ≤ i < j ≤ n where lh ≥ tr, with tr an arbitrary threshold, the structure S is saturated only if no helix H(i, j, lh) with lh ≥ tr can be added to it.

The structure is parametrized because there is a special parameter a that imposes certain restriction on this structure. The space of all saturated secondary structures also gets redefined accordingly. Specifically, from this point on S SAT , resp.

S n SAT contain only saturated structures, resp. saturated structures of size n, that respect the defined parameters. A number of other parameters will be introduced later.

Suppose two helices a ≥ lh and a ≥ lh . Between the helices H(i, j, lh) and H(i , j , lh ), the algorithm assumes two possible relations. First, the helix can be nested in another.

Definition 3.1.4 (Nested helices):

The helix H(i, j, lh) is nested in H(i , j , lh ) if i < i < j < j and there is at least one base w[k] for

w[k] ∈ w[i + a, i -1] ∪ w[i + a, i -1] is unpaired meaning H(i, j, lh
) and H(i , j , lh ) are not the one continuous helix.

Second, the helix can be juxtaposed to another while all of them are nested in the exactly same set of helices.

Definition 3.1.5 (Parallel helices):

The helices H(i, j, lh) and H(i, j, lh) are parallel one to another if either i < j < i < j or i < j < i < j and both helices are, at the same time, either nested in the same helix H(i , j , lh ) or not nested in another helix at all.

The principle of the Saffarian algorithm is, first, to enumerate all possible base pairs. These are then used to construct helices H(i, j, a) that each has the exactly the length of a, the minimum possible helix length. These helices are then assembled to maximum parallel sets.

Definition 3.1.6 (Flat structure):

The flat structure p is the maximum set of parallel helices. Specifically, p(i, j) denotes a maximum set of parallel helices H(i 1 , j 1 , a), . . . , H(i u , j u , a) for a sequence w[i, j], i < i 1 < j 1 < . . . < i u < j u < j. No another un-nested helix H(i , j , a) can be added to it.

Note that by definition, p[i, j] cannot contain any nested elements. For each w[i, j] a certain number of flat structures can be proposed.

Definition 3.1.7 (Flat structures set):

A flat structure set P(i, j) is a set of all possible flat structures p(i, j).

Using the dynamic programming, the flat structures are constructed as follows.

We borrow the expression from the related paper [START_REF] Saffarian | RNA locally optimal secondary structures[END_REF] and update it including the modifications later made by Hélène Touzet [START_REF] Michálik | Efficient approximations of RNA kinetics landscape using non-redundant sampling[END_REF].

• If for i < k ≤ j there is not any helix H(i, k, a), then P(i, j) = P(i + 1, j).

• If such helices exist, then:

P(i, j) =        i<k≤j {H(i, k, a)} ⊕ P(k + 1, j) i<k≤j Sat(H(i, k, a), P(i + 1, j)) (3.1)
The ⊕ operator denotes the concatenation of a base pair to every element of given set, here of the helix H(i, k, a) to every flat structure p of the set P(k + 1, j). The function Sat(H(i, k, a), P(i + 1, j)) checks whether the structure consisting of a concatenation of the helix H(i, k, a) and some p ∈ P(i + 1, j) is saturated. In other words, Sat((i, k), P(i + 1, j)) returns the concatenation of H(i, k, a) and some p ∈ P(i + 1, j) only if there is not another helix that can be added to it.

The next step is to assemble these flat structures p into the saturated secondary structures S SAT . This creates the nested relations between the helices. The assembling step would be done by searching p(1, n) for substring w[1, n] = w, then proceeding with the substrings opened by p(1, n) until such substrings would become too small to host a helix. However, the algorithm proposed as it is is insufficient, since all helices would have length of a or multiples.

To address this issue, if we check w[i, j] that is nested within a H(i -a + 1, j + a -1, a) from some flat structure, we need to suppose the case where an extension of such helix is possible. In that case we suppose an extension by the pair (i, j), elongating H(i -a + 1, j + a -1, a) by one base pair. This makes possible create any helices with lh ≥ a.

It is important to notice that if set P(i, j) of p is nested within the helix H(i -a + 1, j + a -1, a) we have to exclude the helix H(i, j, a) from it. In the opposite case the helices of length of multiples of a could be decomposed in more ways than one, violating the condition of unambiguity.

The algorithm is summed up on Figure 3.1, the substring w[i, j] is specifically treated whether the w[i -1] and w[j + 1] are paired or not.

The problem with this algorithm is its time complexity. While its memory complexity is limited to O(n 2 ), since we need only two-dimensional matrices to store the results, the time complexity depends on the number of flat structures that may be exponential due to the the upper limit for helices being n/2. However, if the number of possible helices for given interval is limited, the complexity becomes polynomial. Therefore, supplementary parametrization can help us to reduce the overall complexity of the Saffarian algorithm.

A i,j

i j i i 1 i 1 + a -1 i 1 + a j 1 -a j 1 -a + 1 j 1 [...] i x i x + a -1 i x + a j x -a j x -a + 1 j x j i i+ a -1 i+ a j-a j-a + 1 j B i,j i j i i 1 i 1 + a -1 i 1 + a j 1 -a j 1 -a + 1 j 1 [...] i x i x + a -1 i x + a j x -a j x -a + 1 j x j i i+1 j-1 j

Parametrization of the Saffarian Algorithm

We mentioned the possibility to parametrize the Saffarian algorithm to search only for saturated secondary structures that satisfy the given conditions and its possible impact on the overall complexity of the algorithm. Two of them, θ and a, were already mentioned. The parametrization of Saffarian algorithm were introduced by Hélène Touzet in a previous modification of the original version of the algorithm.

The list of all parameters is as follows:

• θ θ θminimum length of a base pair. This is due to steric strain and angular tension between the bonds which does not allow for a close nucleotides to be paired.

• a a aminimum length of a helix. The base pair stacks contribute to the stability of the secondary structures, therefore lone base pairs are uncommon [START_REF] Ivo | RNA secondary structures[END_REF].

• π π πmaximum number of unpaired bases in a hairpin. The longer unpaired stretches are less probable due to the condition of saturation.

• m m m is the maximum cumulated size of the unpaired stretches within an internal loop or a bulge. This is the maximum limit of the sum of the lengths α and β of both unpaired stretches. The large bulges/internal loops are considered to be energetically unfavorable [START_REF] Mccaskill | The equilibrium partition function and base pair binding probabilities for RNA secondary structure[END_REF].

• ν ν νmaximum length of a base pair. Specifying this value of this parameter restricts the long-range base pair interactions.

• γ γ γmaximum number of the branches within the multibranch loop. The fixed number of the maximum number of the branches decreases the complexity of the Saffarian algorithm from exponential to polynomial.

All parameters are illustrated on Figure 3. The parameters defined within the Algorithm of Saffarian. ν denotes the maximum base pair span, which allows to define whether we allow the long-rande interactions or not. θ represents the minimum length of a base pair, while π is the maximum length of a hairpin. γ denotes the maximum number of the branches within the multibranch loop.

The advantage of the parametrization of Saffarian algorithm is it does not induce any supplementary complexity cost. As seen on the example of helix length a and minimum base pair length θ, the parameters are directly implemented into the algorithm with and used with the overall complexity O(1). Therefore the space S n SAT can be reduced using the parametrization to restrict the secondary structures without inducing any supplementary complexity cost.

We also redefine S n SAT , resp S SAT as the set of all saturated secondary structures, resp. those with the length n, that fulfill the conditions imposed by parameters a, π, m, ν and γ.

We specifically discuss the minimum helix length a and the maximum number of the branch of the multibranch loop γ. Limiting both of he values can greatly speed up the algorithm of Saffarian at the cost of a loss of the potentially important part of S n SAT . In the case of the minimum helix length a, the default value is 3. While it has been shown that lone base pairs are a rare occurrence [START_REF] Paul | RNA secondary structure: physical and computational aspects[END_REF], excluding the helices of the length of 2 may have bigger impact. For this reason, we analyze the data of experimentally determined secondary structures and study the helices found within them. The similar analysis will be made to the number of the branches within the multibranch loop.

For this purpose, we used the RNA STRAND database [START_REF] Andronescu | RNA STRAND: the RNA secondary structure and statistical analysis database[END_REF] -a database of known experimentally determined secondary structures. We performed its statistical analysis in a following manner:

• Create 5 groups of sequences from RNA STRAND for which the secondary structures, respectively containing all sequences shorter than 140, 300, 400 and 500 nt, the last group containing all sequences.

• For each group, determine a number of helices having at least a certain length.

• For each group, compute the proportion off all helices having least certain length.

The similar analysis is made for the multibranch loop, except we computed the proportion of all secondary structures which multibranch loop with most branches has their given number. The results of both analyses is shown on Figure 3.3A and B respectively.

For the number of the branches within the multibranch loop, we can observe that the structures having having multibranch loops with more than 6 branches are rare. Note that for RNA STRAND, the number of branches counts in also the pair that closes the multibranch loop, while the parameter of Saffarian algorithm does not, so the number of branches in Figure 3.3A corresponds to γ + 1. It is safe to parametrize the Saffarian algorithm with γ = 5, as it leads to a loss of only a negligible proportion of the secondary structures. That would reduce its overall complexity to polynomial, since the maximum number for |p| would be dependent on this parameter instead of the length of the sequence. Such approach seems reasonable considering the number of lost secondary structures is minimal.

The second analysis shows that the helices H(i, j, lh) with the length lh = 1, or lone base pairs (LP), are also rarely observed even among the shorter sequences.

In fact, shorter sequences have lower proportions of LP simply because longer sequences have more base pairs, and therefore also LP. However, around 90% of all helices have at least two or more base pairs. On the other hand, only 75% has H(i, j, lh) with lh ≥ 3, meaning that setting a will lead to a loss of 1/4 of all helices. However, assuming that the functional structures will not be affected too profoundly by exclusion or eventual elongation of these helices, we can allow to set a = 3.

Note that the structures considered here are not all locally optimal or saturated, therefore the statistics for S n L , resp. S n SAT might differ a bit.

From this point on, we suppose that a = 3. The values of γ are unlimited unless noted otherwise.

Computing the Partition Function

The problem of sampling locally optimal secondary structures was already resolved by A. Lorenz and P. Clote [START_REF] William | Computing the partition function for kinetically trapped RNA secondary structures[END_REF] in polynomial time. While the Saffarian al-gorithm samples in exponential time without any parametrization, this exponential parameter is small and therefore the exponential increase is not as remarkable. Plus, because of the parametrization without complexity overhead, the Saffarian algorithm is a valuable tool for searching the locally optimal structures within the defined parametrized space. Since such conditions can be specified, it is easier to look for the structures whose specifics are known. For this reason we decided to work with this algorithm.

Unfortunately, the algorithm of Saffarian does not compute the energies of the secondary structures, instead it enumerates them exhaustively. Therefore, it is impossible to sample the structures since the partition function cannot be computed without energies. Our task is to analyze how to compute the energies of the flat structures so we can compute the partition function using the dynamic programming.

The decomposition of the secondary structures in flat structures allows is to implement the Turner Energy Model by associating the specific flat structures to specific loops. First, consider all possible cases without an extension of an existing helix of some flat structure. We know each such flat structure p contains u helices H(i x , j x , a) such that 1 ≤ x ≤ u. Consequently, their energy can be computed by

E H(i,j,lh) = a-1 y=0
E St (i + y, j -y).

Definition 3.1.8 (Flat structure energy): he energy of the flat structure E p can be computed as

E p = C p + u x=1 lh-1 y=0 E St (i x + y, j x -y) (3.2)
where C p is an energy term that depends on the properties of p.

The properties determining C p are: -number u of helices H(i x , j x , a) such that 1 ≤ x ≤ u, lh x > a; -whether it is nested in another flat structure.

Depending on these properties, we can identify the type of loop that p represents and associate a corresponding value to C p

• If u ∈ N and p is not in another helix, it is an external loop -C f = 1.

• If u = 0 and p is in another helix, p is an hairpin -C p = E H (i, j).

• If u = 1 and p is nested in another helix, p constitutes an internal loop or a bulge depending on i, i 1 , j 1 , j, therefore C p = E ILG (i, i 1 , j, j 1 ).

• If u = 2 and p is nested in another helix, p constitutes a multibranch loop and C p = a + b × u + c × N unp , according to the Equation 2.6.

The extension is a specific case that is equivalent to an elongation of an existing helix, adding an energy term to it.

Definition 3.1.9 (Extension energy): he energy of an extension p E of helix H(i, j, lh) such that H(i, j, a) ∈ p, lh ≥ a is

E p E = E St (i + lh -1, j -lh + 1).
Once each p has been associated to the specific loop case, we can compute the energy of entire secondary structure:

E S = p∈S E p + Ep E ∈S E p E .
With the energies of p, resp. p E being acquired, we can proceed to computation of the Boltzmann distribution. Such a distribution prioritizes the low-energy structures but does not completely ignore high-level energies.

The partition function is computed by employing a DP scheme that recursively includes one flat structure in another and extends their helices. Due to the presence of an extension, we need to distinguish two cases:

• the interval w[i, j] is not included in some helix H(i -lh, j + lh, lh) with lh > a; the partition function is denoted by

Z U i,j ; • the interval w[i, j] is included in a helix H(i -lh, j + lh, lh); the partition function is denoted by Z N i,j .
The values of these two partition functions can be computed as:

Z U i,j = p(i,j)∈P(i,j) e -E p (i,j) k B T × x c=1 Z H ic+a,jc-a Z H i,j = p(i,j)∈P(i,j)/H(i,j,a) e -E p (i,j) k B T × x c=1 Z H ic+a,jc-a + e -E p E (i,j) k B T × Z H i+1,j-1 (3.3)
In the case of Z B i,j , we need to exclude H(i, j, a) since it would extend the helix H(i -a, j -a, a) of the parent flat structure p and would create ambiguous situation since such decomposition is handled by the extensions. The extension is then included in Z H i,j since we know that it is included in some H(i -a, j -a, a). The partition function Z is given by Z H 1,n .

The computation of the energies of the loops in RNANR employing Saffarian algorithm uses the functions from VIENNARNA which also includes all necessary energy parameters. These functions as well as DP scheme from Equation 3.3 are implemented into the RNANR in C language. The stochastic backtrack constituting the base of the sampling in the way similar to Section 2.5.6.

The similar approach was used independently in works of Jérôme Waldispühl and Peter Clote [START_REF] Waldispühl | Computing the partition function and sampling for saturated secondary structures of RNA, with respect to the Turner energy model[END_REF].

Formalization

Here we express the DP scheme employed in Saffarian algorithm using the formalism we established in Section 2.6. Its formalization is quite simple. First, we need to define states Q. As with Zuker and McCaskill algorithm, w[i, j] alone is not sufficient to define a state, but we need to know whether (i -1, j + 1) ∈ S SAT (i, j), with S SAT (i, j). a saturated structure constructed by DP scheme. We therefore define two types of states:

• q U i,j = w[i, j]|(i, j) / ∈ S SAT (i, j) • q H i,j = w[i, j]|(i, j) ∈ S SAT (i, j)
The set of all states is then defined by

Q Saf = {q U i,j , q H i,j | 1 ≤ i < j ≤ n}.
There are two types of derivation and associated constructors:

• λ flat i,j,p(i,j) -associates a flat structure p(i, j) ∈ P(i, j) to w[i, j] • λ ext i,j -an extension of a helix H(i -lh, j + lh, lh), lh > a

We can consequently define ρ Saf :

ρ Saf (q U i,j ) = p(i,j)∈P(i,j)
({q H i 1 +a,i 1 -a , . . . , q H iu+a,iu-a }, λ flat i,j,p(i,j) , E p (i, j))

ρ Saf (q H i,j ) =                p(i,j)∈P(i,j)/ /H(i,j,a)
({q H i 1 +a,i 1 -a , .., q H iu+a,iu-a }, λ flat i,j,p(i,j) , E p (i, j))

({q H i,j }, λ ext i,j , -E p E (i, j)) (3.4)
Here q root = q U 1,n . We have thusly completely defined the DP Scheme used in Saffarian algorithm as (Q Saf , q U 1,n , ρ Saf ).

To compute the partition function, we proceed by application of the Equation 2.6.2. This gives us:

PfS Saf (q U i,j ) = q U i,j λ flat i,j,p(i,j) -----→ E p (i,j) {q H i 1 ,j 1 ,...,q H iu,ju } e E p (i,j) k B T × u i=1 PfS Saf (q H ix,jx ) PfS Saf (q H i,j ) = q H i,j λ flat i,j,p(i,j) -----→ E p (i,j) {q H i 1 ,j 1 ,...,q H iu,ju } e -E p (i,j) k B T × u i=1 PfS Saf (q H ix,jx ) + + e E p (i,j) k B T × PfS Saf (q H i+1,j+1 )
Finally, the search space for each state q can be defined as:

L Saf (q U i,j ) = q U i,j λ flat i,j,p(i,j) -----→ E p (i,j) {q H i 1 ,j 1
,...,q H iu,ju } λ p (L(q H i 1 +a,i 1 -a ), . . . , L(q H iu+a,iu-a ))

L Saf (q H i,j ) =                q U i,j λ flat i,j,p(i,j) -----→ E p (i,j) {q H i 1 ,j 1
,...,q H iu,ju } λ p (L(q H i 1 +a,i 1 -a ), . . . , L(q H iu+a,iu-a ))

λ ext (L Saf (q H i+1,j-1 )) (3.5)
The stochastic backtrack is then formalized by using the Equation 2.18, where f λ is either the function computing the partition function for applying p or an extension of one of helices of p .

Comparison of Saffarian and Turner Local Minima

Saffarian algorithm assumes S SAT is locally optimal if it is saturated. However, nothing so far confirms such assumption, specifically with regards to Turner Energy Model. Since we use this energy energy model to compute the energies of all flat structures and secondary structures determined by the Saffarian algorithm, we would like to study whether the saturated secondary structures are effectively the local minima according to it.

The local optimality implies that all neighbors S neigh ∈ V S SAT to S L satisfy:

E S SAT < min Sx∈V S SAT E Sx .
To validate this condition we apply the method of gradient walk. This method searches for all possible neighbors of a saturated secondary structure S SAT and calculates their free energy. Should one of neighbors have lower energy than S SAT , then the neighbor becomes then new minimum and the entire process is repeated until we reach the true local minimum. If S SAT is local minimum then no descent will be possible.

For this purpose, we employ the gradient walk function implemented in VIEN-NARNA.

It must be pointed out that the space of secondary structures accessible by the gradient walk is not the same as the one accessible by the Saffarian algorithm. While the gradient walk can access the entirety of S n , the Saffarian algorithm can access only the parametrized space S n SAT ⊂ S n . This means that after performing the gradient walk on S SAT we might find ourselves with S / ∈ S n SAT . It is important to know whether the structure is not locally optimal because S n SAT does not include such structure or it is caused by the conception of the algorithm. For this reason here we will separate the results obtained for the gradient walk to those that are within search space, or where the final S is in parametrized S n SAT , and those that are not not, or are outside search space.

To compute these values, we used 154 sequences from RNA STRAND, or all sequences from the database version from 04/11/2011 having between and including 120 and 170 nt. This ensures the sufficient variability for the secondary structures and avoiding the situations where saturated secondary structures would be also locally optimal by random chance.

The proposed pipeline of validation of saturated structures as locally optimal is:

1. Using the Saffarian algorithm, sample 1000 unique saturated secondary structures S SAT . This can be either done by sampling S SAT and removing duplicates until we have 1000 uniques, or more easily performing the nonredundant sampling that will be detailed in the next section.

2. For each S SAT , perform a gradient walk using a dedicated VIENNARNA function.

3. Compare the gradient walk result S grad with corresponding S SAT .

4. For the S grad that are different, check out whether they satisfy the parametrization of Saffarian Algorithm, ie. they validate the parameters θ, a, π, m, ν and γ.

The results are listed in Table 3.1. The S grad that are identical to starting S SAT are counted in the category Within search space. We found that of 59.57% of S grad within search space, of which about 90% are locally optimal according to the Turner Energy Model. For the S grad ∈ S SAT that differ from S SAT , the true local minimum is within two base pairs most of the time (average at 1.550), and within the 2 kcal.mol -1 of free energy (average at ∆∆1.258). The saturated secondary structures S SAT are therefore identical to S grad in 52%, and in most of others S grad is quite close. However, in the majority of cases where S grad and S SAT differ we have S grad / ∈ S SAT . This suggests that the local optimum cannot be reached because it is outside S SAT . On the other hand, there are some case where the saturated structure is not locally optimal. However, since such results are few and far apart, we can indeed assume that S SAT ∈ S L and therefore that Saffarian algorithm returns locally optimal structures. These structures are identical to S SAT in about 50% of cases. The category denoted 'Outside search space' specifies locally optimal structures obtained the same way and not fulfilling the parametrization of the Saffarian algorithm. ∆∆G specifies the average difference of energy between the two structures, while the average base pair distance accounts to the average number of different base pairs between them.

Samples%

Note that this experiment verifies that S SAT are locally optimal with respect to the Turner Energy Model, but it does not confirm whether the space of saturated secondary structures completely covers the space of locally optimal structures that respect the given parameters. This is not problematic though, since in this case we are interested only in saturated secondary structures.

The Non-redundant Sampling Problem 3.2.1:

• Input: DP scheme that can be formalized as (Q, q root , ρ);

• Output: set of length N of unique secondary structure samples in normalized Boltzmann distribution.

The advised non-redundant sampling strategy has to fulfill two conditions:

• The method must be applicable to any DP scheme computing the partition function for the purpose of predicting suboptimal secondary structures of RNA.

• The method should be easy to implement, therefore constituting an interface to the adjacent DP scheme with the lowest number of interactions possible.

The general concept of non-redundant sampling was introduced and discussed for weighted generation from context-free grammars by William Lorenz and Yann Ponty [START_REF] Andrew | Non-redundant random generation algorithms for weighted context-free grammars[END_REF]. We push said method further, adapt the result for the purpose of RNA secondary structure prediction and put it into practice.

Here, we will, after introducing the necessary notions and explaining the general principle of non-redundant sampling, research the method that fulfills aforementioned criteria. We will discuss on the points that must be addressed to create such method. The two main interdependent issues is the memorization process of the sampled secondary structures and their decomposition and how to ensure that the distribution of non-sampled structures stays unbiased. While the first experiments will be made on the RNANR, we will try to generalize the principle as much as possible using the formalism we introduced except for specific examples.

We will follow by the experiments demonstrating the efficiency of this approach when compared to RNA usual sampling methods. We will then describe the implementation of researched non-redundant sampling method in the most optimal way and backing the general applicability of it to the current DP schemes used for sampling of secondary structures from Boltzmann distribution.

Parse tree and Incomplete structures

Before explaining the non-redundant sampling, we must introduce the notion of parse tree and incomplete structures in the context of secondary structure sampling. Most DP schemes only consider what happens for the given q at the given point, but not the events taking place outside of it. The principle of nonredundant sampling necessitates to know this information, and to express it from that generated by DP scheme. The specified concepts are introduced for this reason.

Each DP scheme (Q, q root , ρ) can be captured as a context-free grammar [START_REF] Andrew | Non-redundant random generation algorithms for weighted context-free grammars[END_REF], where

• The alphabet is the set D of constructors λ;

• The non terminal symbols are the states q ∈ Q;

• The production rules are derivations q λ -→ c λ {q 1 , . . . , q u }, taking q and returning λ and q 1 , . . . , q u -for now we ignore the score contribution c λ ;

• The axiom is q root . Therefore each secondary structure can be represented as a parse tree.

Definition 3.2.1 (Parse tree):

A secondary structure parse tree is a tree T(w, Q, q root , ρ) for which we have:

• Internal nodes λ ∈ D -each parent of u u u subtrees T 1 , . . . , T u with derivation q λ -→ c λ {q 1 , . . . , q u }; • Leaves λ term ∈ D such that q λ term ----→ c λ term ∅.
Such a tree represents a structure where an application of constructor leads to application of λ 1a and λ 1b etc., each building a specific element. Now suppose the structure is not completely derived yet on states q. In that case some of the leaves will not be represented by a constructor λ, but instead by some state q which can be derived by a number of derivations q λ -→ c λ {q 1 , . . . , q u }.

Definition 3.2.2 (Immature parse tree): Let Q N be the list of underived states. An immature secondary structure parse tree is a tree T I (w, Q, Q N , q root , ρ) for which we have:

• Internal nodes λ ∈ D -each parent of u u u subtrees T 1 , . . . , T u with derivation q λ -→ c λ {q 1 , . . . , q u };

• Leaves λ term ∈ D such that q λ term ----→ c λ term ∅ or q ∈ Q N .
The immature parse tree evolves as the states q get derived as λ term ∈ D such that q

λ term ----→ c λ term ∅.
We can represent this evolution by a decomposition tree.

Definition 3.2.3 (Decomposition tree):

A decomposition tree is a tree T (w, Q, q root , ρ)

for which:

• Each internal node V is a immature secondary structure parse tree T I (w, Q, Q N , q root , ρ) .

• Each leaf V term is mature secondary structure parse tree T(w, Q, q root , ρ).

• Each edge V 1 → V 2 represents a derivation λ term ∈ D such that q λ term ----→ c λ term
∅, where q is a leaf of a parse tree in some node V.

Example 3.2.1. Suppose a sequence w and a DP scheme (Q Ex , q root , ρ Ex ) where Q ex = {q root , q 1A1 , q 1B1 , q 1B1 , q 1B2 } and the ρ Ex is given by

ρ ex (q root ) = {({q 1A1 , q 1A2 }, λ 1A , c 1A ), ({q 1B1 , q 1B2 }, λ 1B , c 1B )} ρ ex (q 1A1 ) = {(∅, λ 2A , c 2A )} ρ ex (q 1A2 ) = {(∅, λ 2B , c 2B )} ρ ex (q 1B1 ) = {(∅, λ 2C , c 2C )} ρ ex (q 1B2 ) = {(∅, λ 2D , c 2D )}
The decomposition tree is given by Figure 3.4. Each internal node is represented by an immature parse tree, and each leaf is a mature parse tree. At the beginning there is only the initial state dproot. This state gets derived by one of two available derivations,

q root λ 1A ---→ c λ 1A {q 1A1 , q 1A2 } and q root λ 1B ---→ c λ 1B
{q 1B1 , q 1B2 }, giving two new different states for each case. We derive the states q 1A1 and q 1A2 that follow

q root λ 1A ---→ c λ 1A
{q 1A1 , q 1A2 }, resp. q 1B1 and q 1B2 that follow

q root λ 1B ---→ c λ 1B
{q 1B1 , q 1B2 } consecutively, leading to the mature parse trees at the end of each branch.

The important point is the order in which the states q of an incomplete structure are derived. This order must be deterministic, otherwise multiple decomposition trees may be created for the same DP scheme and sequence. In our case we assume the state q to be derived is the rightmost one, and that all trees for the same DP scheme (Q, q root , ρ) and sequence w have always the same layout.

q r oot λ 1A q 1A1 q 1A2 λ 1A λ 2A q 1A2 λ 1A λ 2A λ 2B λ 1B q 1B1 q 1B2 λ 1B λ 2C q 1B2 λ 1A λ 2C λ 2D Figure 3
.4: Decomposition tree T (w, Q ex , q root , ρ ex ). Each internal node is constituted from an immature parse tree T I (w, Q ex , Q N , q root , ρ ex ), while the leaves represent a mature parse trees T I (w, Q ex , q root , ρ ex ) and T(w, Q ex , q root , ρ ex ). The nodes are connected by the derivations q λ --→ c λ q 1 , . . . , q u on the nodes q of immature trees.

We compute the search space L(T) associated to a tree T(w, Q, q root , ρ) or T I (w, Q, q root , ρ) in a recursive manner:

L(T) =    λ(L(T 1 ), • • • , L(T u )))
if λ an internal node or mature leaf λ(L(q)) if q an immature leaf (3.6) Theorem 3.2.1 (Search space inclusion): Suppose a decomposition tree T where V contains two parse trees T I and T I connected by an edge representing a derivation q λ -→ c λ q 1 , . . . , q u from T I to T I , q ∈ T I . In that case

L(T I ) ⊆ L(T I ).
Proof 3.2.1. The derivation q λ -→ c λ {q 1 , . . . , q u } implies a tree T I contain a state q that can be derived. Since Definition 2.6.1 defines constructor λ as a function with properties of cartesian product and the T I , resp. T I can be defined recursively by the properties of subtrees T 1 , • • • , T u , we can separate L(q) from other terms in L(T I ). The derivation q λ -→ c λ {q 1 , . . . , q u } introduces λ and the space of structures λ(L(q 1 , . . . , q u )). Equation 2.6.1 implies that λ(L(q 1 , . . . , q u )) ⊆ L(q).

Due to the properties of λ as cartesian product, it follows L(T I ) ⊆ L(T I ).

We define an incomplete secondary structure, whose notion is directly related to that of an immature secondary structure parse tree.

Definition 3.2.4 (Incomplete structure):

We call an incomplete structure S I a secondary structure that contains at least one state q that was not yet derived. For a sequence w, an incomplete structure is always of size n.

Example 3.2.2. The illustrated example for DP scheme

(Q Zuk , q F 1,12 , ρ Zuk )
employed in Zuker algorithm is available on

• q C 4,8 λ H 4,8 ---→ c λ H 4,8
∅ that creates a hairpin.

An incomplete structure delimits a certain subspace of S n .

Definition 3.2.5 (Incomplete structure space):

The incomplete structure space S I of S I is a space of secondary structures that can be obtained by deriving on underived state q by q λ -→ c λ {q 1 , . . . , q u }.

Example 3.2.3. Assume that the state q C 4,8 is the only underived state of of an incomplete structure S Ia from Figure 3.5. This state can be derived by q C The partition function of an incomplete structure S I is

S I a = {S 1 , S 2 }. A C C C A A A G G A G A C A A A G C A A A G 1
Z S I = S∈S I Z S .
In the decomposition tree T , each node in V can be seen as an incomplete secondary structure S I where leaves q are underived states in secstrinc. Each node is an immature parse tree T I with known set of underived states q and a chain of performed derivations, as well as a subset of structures L(T I ). Therefore T I represents a secondary structure S I and L(T I ) = S I . This also implies that for two incomplete structures S I a and S I b connected by an edge in T , where S I b was obtained by deriving q of S I a , we have

S I b ⊆ S I a .
A mature tree T(w, Q, q root , ρ) represents a completed structure, since it does not have any leaf q. Said structure is constructed by following a path in T from root to specified leaf.

General Principle of Non-redundant sampling

Suppose two incomplete structures S Ia , represented by an immature parse tree T Ia with a leaf q the next state to be derived, and S Ib , represented by T Ib such that S Ib is obtained by derivation q λ -→ c λ {q 1 , . . . , q u } of q in S Ia and is consequently connected by an edge in an associated decomposition tree T (w, q, q root , ρ). In addition to the Equation 2.21, the probability p(λ | q) can be expressed as

p(λ | q) = Z S Ib Z S Ia (3.7)
If we want to perform a non-redundant sampling, we need to make sure that at every moment we choose a derivation, the newly generated incomplete structure S Ib will not lead to S I composed only from the structures that were already sampled. This must be verified every time a new λ is selected. This is simple if one knows what structures S F were already sampled and which S F ∈ S I , with S F the corresponding forbidden structure. We denote the space of forbidden structures, structures than must not be generated, by F.

In the context of non-redundant sampling, F contains all already sampled sampled structures. The set of forbidden structures F S I that can be generated by completing S I is

F S I = F ∩ S I .
Its partition function Z F S I is

Z F S I = S∈F S I Z S .

Definition 3.2.6 (Non-redundant emission probability):

Let F S I a and F S I b be the forbidden structures that can be generated by completing S Ia , resp. S Ib . The probability p(λ | q, F S I a , F S I b ) of choosing a derivation

q λ -→ c λ {q 1 , . . . , q u } is p(λ | q, F S I a , F S I b ) = Z S Ib -Z F S I b Z S Ia -Z F S I a . (3.8) Theorem 3.

(Non-redundant generation probability):

Using the nonredundant emission probability for each derivation, we generate a structure S ∈ L(q root ) -F with global probability

p(S | F) = Z S -1 S∈F × Z S Z -Z F . (3.9)
where 1 S∈F is an indicator function that takes value 1 if S ∈ F and 0 otherwise. This is effectively the expression of a probability to sample S; if S ∈ F we have p(S | F) = 0, otherwise it scales accordingly by the Boltzmann factors of the structures not in F.

Proof 3.2.2.

From the decomposition tree T (w, q, q root , ρ), the generation of a structure S can be represented by a path of T where leaf T Iv represents S:

T I0 → T I1 → • • • → T Iv = T.
Each T Iv represents an incomplete secondary structure, therefore

S I0 → S I1 → • • • → S Iv = S.
Assume that F S Ix is the set of forbidden structures which can be generated from S Ix where 1 < x < v. The probability p(S | F) of choosing S is given by the probability of selecting each λ for which q λ -→ c λ {q 1 , . . . , q u } leads to it:

p(S | F) = p(λ 1 | S I0 , F S I 0 , F S I 1 ) × p(λ 2 | S I1 , F S I 1 , F S I 2 ) × • • • • • • × p(λ v | S Iv-1 , F S I v-1 , F S I v )
From there, we get

p(S | F) = Z S I1 -Z F S I 1 Z S I0 -Z F S I 0 × Z S I2 -Z F S I 2 Z S I1 -Z F S I 1 × • • • × Z S Iv -Z F S I v Z S Iv-1 -Z F S I v-1 p(S | F) = Z S Iv -Z F S I v Z -Z F p(S | F) = Z S -1 S∈F × Z S Z -Z F . (3.10)
Since this is exactly the expression from Equation 3.9, we validate the Theorem 3.2.1.

The expression is in principle similar to the one employed in Algorithm 1 from the paper presented by Lorenz et al [START_REF] Andrew | Non-redundant random generation algorithms for weighted context-free grammars[END_REF]. Notice that if S Ib = F S Ib , then correctly p(λ | q, F S Ia , F S Ib ) = 0.

Unfortunately, this expression cannot be applied directly. As already pointed out, the DP schemes employed for the RNA secondary structure prediction do not consider what is happening outside their specific state. Thus, they do not compute Z S I but only PfS(q) associated to a certain state q. Therefore we need to research a way how to use the values of PfS(q) to access Z S I to perform a non-redundant sampling. The second problem is that all these values as well as those of F S I need to be stored after the computation and must be easily accessible at given moment. Here we separately address both problems and discuss the solutions.

Computing the probability with forbidden structures

This part discusses how to rely the values from the Equation 3.8 and PfS(q) that is computed using the DP scheme (Q, q root , ρ). We assume tree the decomposition tree T (w, Q, q root , ρ) and an incomplete structure S I , represented by an immature parsing tree T I (w, Q, Q N , q root , ρ) with the following properties:

• Derivations q λ 1 , . . . , q λu , attributing free energies c λ 1 , . . . , c λu , were applied beforehand to their respective states.

• q is a leaf of T I for which we have a set of derivations

q λ -→ c λ {q 1 , . . . , q z }
and according to T (w, Q, q root , ρ) and this state is derived next.

• T I may or may not have other leaves q y , 1 ≤ y ≤ v for which we have derivations q λy ---→ -c λy {q 1y , . . . , q vy } and their derivations follow after that of q .

The partition function of S I can be separated to the contributions c λ of different λ and the contributions PfS(q) of all states to which no λ is attributed.

Theorem 3.2.3 (Incomplete structure partition function):

The partition function Z S I of an incomplete secondary structure S I defined above can be computed by

Z S I = u x=1 e -c λx k B T × PfS(q) × v y=1
PfS(q y ).

(3.11)

Proof 3.2.3. We prove the above by induction on the number of the derivations M necessary to obtain Z S I with underived states q.

Base case (M = 0): If no derivation was performed then the only underived state present is the starting state q root . Therefore:

Z S I0 = PfS(q root )
In this case Z S I0 = Z, since q root is the starting state. We have already shown in Proof 2.6.2 that under the conditions of completeness and unambiguity we have

PfS(q) = s∈L(q) e -c(s) k B T
and since completeness implies L(q root ) = S n , Z S I0 = Z Sn = PfS(q root ).

Induction (M = m): Suppose S I was obtained after m -1 derivations q λx -→ λx {q 1x , . . . , q ux } with 1 ≤ x ≤ m -1. We suppose by induction that the expression is valid for every S I on which m -1 was performed and that the next incomplete structure S Inext , represented by an immature parse tree T I (w, Q, Q N , q root , ρ), is generated by deriving S I with q λ -→ c λ {q 1 , . . . , q z }. This implies S Inext is accessed after m derivations. From Equation 3.11, we have

Z S Inext = e -c λ k B T × u x=1 e -c λx k B T × z a=1 PfS(q a ) × v y=1
PfS(q y ).

(3.12)

From the Equation 2.21, we can express

e -c λ k B T × z a=1
PfS(q a ) = p(λ | q) × PfS(q). Therefore the Equation 3.12 gives:

Z S Inext = p(λ | q) × PfS(q) × u x=1 e -c λx k B T × v y=1
PfS(q y )

Z S Inext = p(λ | q) × Z S I
which is exactly the same expression given by Equation 3.7, validating the induction step.

Consequently, we can express the Equation 3.8 as

p(λ | q, F S I , F S I next ) = e -c λ k B T × u x=1 e -c λx k B T × z a=1 PfS(q a ) × v y=1 PfS(q y ) -Z F S I next u x=1 e -c λx k B T × PfS(q) × v y=1 PfS(q y ) -Z F S I
.

This expression cannot be easily simplified. However, we can reduce it by Z S I / PfS(q):

p(λ | q, F S I , F S I next ) = PfS(q) PfS(q) × e -c λ k B T × u x=1 e -c λx k B T × z a=1 PfS(q a ) × v y=1 PfS(q y ) -Z F S I next u x=1 e -c λx k B T × PfS(q) × v y=1 PfS(q y ) -Z F S I p(λ | q, F S I , F S I next ) = Z S I PfS(q) × e -c λ k B T × z a=1 PfS(q a ) -Z F S I next Z S I -Z F S I p(λ | q, F S I , F S I next ) = e -c λ k B T × z a=1
PfS(q a ) -PfS(q)

Z S I × Z F S I next Z S I -PfS(q) Z S I × Z F S I (3.13)
The final equation gives us the value that is similar to the general stochastic expression from Equation 2.21. Since the proportions between the structures outside F are conserved, the probability p(λ | q, F S I , F S I next ) is unbiased. This is the expression that allows by which we compute the probability of choosing λ in a non-redundant manner.

The values of Z F S I can be accessed from a specific data structure that stores them due to the necessity to know which subset of F can be generated at given moment. This structure will be detailed later. We still need to compute Z S I as well. This can be done recursively, as we know that

Z S I0 = Z and Z S Inext = Z S I × p(λ|q)
and p(λ|q) without forbidden terms can be computed from Equation 2.21. The exact principle can be,when using the decomposition tree T (w, Q, q root , ρ) with sequence w and the DP scheme (Q, q root , ρ), summed up in points as follows:

1. Initialize Z S I → Z and q → q root .

2. Compute p(λ | q root , F S I , F S I next ) for every λ such that q root λ -→ c λ {q 1 , . . . , q u } 3. Randomly choose q root λ -→ c λ {q 1 , . . . , q u } depending on the probabilities.

4. Compute p(λ|S I ), then compute Z S Inext .

5. Set the value of Z S I to Z S Inext .

6. Set q root = q 1 and store other states in T (w, Q, q root , ρ). Repeat steps 2 -5 until S I = {S}.

Once S is completed, it is added to F. Such approach demands the minimum changes to the original DP scheme, since only updating Z secstrinc and subtracting the values of F S I is necessary. Both require a specific data-structure which is discussed next.

Adjacent data-structure for non-redundant sampling Problem 3.2.2:

• Input: Decomposition tree T (w, Q, q root , ρ).

• Output: a data-structure d that:

-Adds nodes that stores for each T I (w, Q, Q N , q root , ρ) representing S I a value Z forb S I if the node does not exist and S I was generated.

-Accesses Z forb S I if a new structure S ∈ S I is generated and updates Z forb S I at the end of given generation cycle by adding Z S to it.

The problem implies the adjacent data structure d must be based on T . Storing the T (w, Q, q root , ρ) with all nodes would be a memory consuming task. However, each node T I of T (w, Q, q root , ρ) can be reduced only to the derivation that was performed last at the given point (Figure 3.6) or an axiom, dproot. We therefore need to store, along with Z forb S I , only λ.

Definition 3.2.7 (Decomposition tree reduction): Let T be the space of all trees T I (w, Q, Q N , q root , ρ). The decomposition tree reduction is a function

R : T → D
that associates to each node T I of T a λ such that q λ -→ c λ {q 1 , . . . , q u } was the last derivation performed in T I .

q r oot λ 1A q 1A1 q 1A2 λ 1A λ 2A q 1A2 λ 1A λ 2A λ 2B λ 1B q 1B1 q 1B2 λ 1B λ 2C q 1B2 λ 1A λ 2C λ 2D q r oot λ 1B λ 2C λ 2D λ 1A λ 2A λ 2B Figure 3
.6: Simplification of a decomposition tree T (w, Q, q root , ρ) into a datastructure d(w, Q, q root , ρ). Each node is reduced to the derivation that was performed last.

Definition 3.2.8 (Non-redundant sampling data structure):

The datastructure for non redundant sampling is a tree d(w, Q, q root , ρ), based on a decomposition tree

T (w, Q, q root , ρ) = G(V, E)
where

• E is the set of edges from T ; • V is a set of nodes R(T I ) with T I ∈ T .
Constructing the entirety of d is memory-consuming even with the application of R. However, d can be constructed in parallel with the non-redundant sampling -if a node λ = R(T I ) that would represent T, resp. S I is absent, it means no structure S ∈ S I was sampled yet, therefore Z F S I = 0. If the node λ does exist, the value Z F S I it stores is used to compute p(λ | q, F S I , F S I next ) via the Equation 3.13.

The probabilities are then used to sample which λ will be attributed to q, and if the node does not exist yet it is created with y∈F S I Z y = 0.

Note that the reduction R is accompanied by the loss of the information on the layout of T and specifically the order of derivation of states q since there may more of them at once. However, the application of recursive function handles this problem as well as using a stack W storing the nodes that must be treated and popping them one at a time.

Once we complete a structure and generate S, equivalent to obtaining a mature parse tree T, we traceback the path from leaf to root while adding Z S to the value Z F S I of each node λ = R(T I ). We have shown that in T every structure of children node T I is represented by its parent T I , therefore S is contained by all nodes from its respective mature parse tree T to root.

A special mention goes to the linking of the children and parent node. This is important and has the impact on the efficiency of the non-redundant sampling implementation. The obvious method would be to create a hash table where the values of λ identifying each node serve to compute the hashing function. However, the simple linked list seems to be more efficient and is simpler to implement.

In that case the nodes are liked in the same order as they are treated in DP scheme (Q, q root , ρ). In that case we move to the next node of linked list it corresponds to the case in DP scheme, and if it does not, it means that no structure from S I represented by the node labeled λ was treated beforehand, implying Z F S I = 0. If such a structure is sampled, the node is inserted in linked list at the corresponding position, preserving the order.

Example 3.2.4. The concrete example is given at Figure 3.7 for DP scheme

(Q Zuk , q F 1,9 , ρ Zuk )
and sequence w. Here no structure containing a pair (2, 8) was sampled beforehand, meaning node λ pair 1,8,2 was not created before. If such structure is to be sampled, respectively λ pair 1,8,2 is selected with p(λ pair 1,8,2 | q F 1,8 , F S I 1 , ∅) with S I 1 being the incomplete structures represented by node λ pair 1,9 , the node is inserted at the position copying the order in DP scheme.

λ unp 1,9 ZF SI 1 q F 1,8 q F 1,8 q F 1,8 ({q F 1,7 }, λ unp 1,8 , 1) q F 1,7 q F 1,7 q F 1,7 q C 4,7 q C 4,7 q C 4,7 ({q F 1,2, q C 4,7}, λ pair 1,8,3 , EBP 3,8 ) λ unp 1,8 ZF SI 21 λ pair 1,8,3 ZF SI 23 q F 1,8 q F 1,8 q F 1,8 λ unp 1,9
ZF SI 1

({q F 1,7 }, λ unp 1,8 , 1)
({q F 1,1 , q C 3,7 }, λ pair 1,8,2 , EBP 2,8 ) ({q F 1,2 , q C 4,7 }, λ pair 1,8,3 , EBP 3,8 ) λ unp 1,8 ZF SI 21 λ pair 1,8,2 ZF SI 22 = 0 λ pair 1,8,3 ZF SI 23 q F 1,7 q F 1,7 q F 1,7 q C 3,7 q C 3,7 q C 3,7 q C 4,7 q C 4,7 q C 4,7 A C C A A A G G A A C C A A A G G A A C C A A A G G A A C C A A A G G A q F 1,8 q F 1,8 q F 1,8 q F 1,7 q F 1,7 q F 1,7 q C 3,7 q C 3,7 q C 3,7 q C 4,7 q C 4,7 q C 4,7 ({q F 1,7 }, λ unp 1,8 , 1) ({q F 1,1 , q C 3,7 }, λ pair 1,8,2 , EBP 2,8 ) ({q F 1,2 , q C 4,7 }, λ pair 1,8,3 , EBP 3,8 )
Figure 3.7: Linking the parent and children nodes in the data-structure for nonredundant sampling in tree d(w, Q Zuk , q F 1,9 , ρ Zuk ). Here no secondary structure with pair (2, 8) was created beforehand, therefore node λ pair 1,8,2 is not yet created. If such structure is sampled the not gets inserted to its respective position, copying the order within DP scheme (Q Zuk , q F 1,9 , ρ Zuk ).

Non-redundant sampling algorithm

The algorithm of non-redundant sampling algorithm for DP scheme (Q, q root , ρ) consists of the data structure d(w, Q, q root , ρ) storing the informations to compute the probabilities of choosing non-redundantly a certain case for the state q according to the relations introduced in Section 3.2.2. Here we resume the entire algorithm. These probabilities evolve with each sampling iteration implying the necessity of update of informations. The d itself is constructed as different derivations are performed. Its pseudocode can be consulted on Algorithm 31.

Assume the number of the unique samples we want to generate is N. The principle of the non-redundant sampling can be resumed as follows:

1. Pre-compute all DP matrices of DP scheme (Q, q root , ρ).

2. Initialize the adjacent structure for non-redundant sampling d(w, Q, q root , ρ)

by placing q root at its root.

3. Set Node = q root , Z S I = Z and prepare a empty stack W.

For each sample:

(a) Get the value of Z F S I q stored in the Node.

(b) List all derivations q n λ --→ c λ {q 1 , . . . , q w } where q n is:

• q root if Node = q root ;
• first node from stack W if w = 0;

• q 1 , given by q λ -→ c λ {q 1 , . . . , q u } otherwise.

(c) Initialize a sum Z temp = 0.

(d) Check if node λ is present. Extract Z F S I from it if is, otherwise Z F S I = 0. Compute p(λ a | q n , F S I , F S I ) (Equation 3.13).
(e) Update Z temp = Z temp + p(λ | q, F S I q , F S I ). If Z temp > r, choose derivation q λ --→ c λ {q 1 , . . . , q w } and update Z S I = Z S I ×p(λ | q, F S I , F S I ).

(f) If for q λ --→ c λ {q 1 , . . . , q w } we have w > 1, put all states other than q 1 in W.

(g) Create node λ linked to Node with Z F S I q 1 = 0 if it does not exist. Set Node = λ .

(h) Repeat steps (a) -(g) until no derivation is possible.

(i) Backtrack from q to q root while updating F S I q for all nodes λ on path. Repeat substeps of step 4 for each sample.

The time complexity of this algorithm is that of the adjacent DP scheme (Q, q root , ρ), since for each decision in DP scheme it performs a limited number of steps. Its memory complexity is O(|Q| + |D|) due to necessity to remember each choice of derivation that was performed.

100

Notice that the main interaction between the DP scheme (Q, q root , ρ) and the layer assuring the non-redundancy modeled by d(w, Q, q root , ρ) is only when extracting F S I , creating new node λ if absent and backtracking once the complete structure is generated. This makes it very easy to implement it into DP scheme datastructure. However, the non-redundant layer must be adapted for the DP scheme in question.

Function Main():

Data: A DP scheme (Q, qroot, ρ),

RNANR: Implementation and results

The non-redundant sampling was implemented into the the program using the Saffarian algorithm and sampling the locally optimal secondary structures -RNANR, standing for "RNA non-redundantly". The general principle stays the same, the only specificity is the flat structures p. Besides the non-redundant sampling and the computation of probabilities using Turner Energy Model researched and implemented by us, it also includes the exhaustive enumeration and parametrization developped by A. Saffarian et al [START_REF] Saffarian | RNA locally optimal secondary structures[END_REF] and later modified by H. Touzet [START_REF] Michálik | Efficient approximations of RNA kinetics landscape using non-redundant sampling[END_REF].

The extension by non-redundant sampling was coded like the original program in C. The most energy-computing functions related to the Turner Energy Model are computed by using VIENNARNA library.

The RNANR software is available at the site: https://project.inria.fr/ rnalands/software/rnanr/.

This section discusses the various aspects of the software, the problems that were encountered, and the results of the different tests that were realized.

Numerical instability issues

The problem related specifically to the non-redundant sampling is the numerical instability. This is not the problem for the regular sampling since the contribution of the sampled structures are not cumulated, but in the case of the non-redundant sampling, these values are summed up and if the sampling performed is of the sufficient size, the Boltzmann factors Z S of the secondary structures S become eventually too small for the small precision deviations to become an issue. This leads to selecting wrong or malformed secondary structures once Z S reaches values that are numerically indistinguishable from zero. This is more probable to happen for bigger number of unique samples, and is especially noticeable for the sequences having secondary structures important free energy differences -specifically for the the artificially designed riboswitch SV11 generated by Qβ replicase [START_REF] Biebricher | Structural analysis of self-replicating RNA synthesized by Qβ replicase[END_REF] the numerical instability appears as soon as after 5000 samples.

To counter the described problem, we employ the GNU MPFR library [START_REF] Fousse | MPFR: a multiple-precision binary floating-point library with correct rounding[END_REF], a multiple-precision floating-point computation library based on GMP [START_REF] Granlund | GNU MP: The GNU Multiple Precision Arithmetic Library[END_REF] that allows to store floating point values with arbitrary precision. This includes the values Z S that are too small when compared to Z, for which higher precision is used. In the case of RNANR, the values of partition functions and Boltzmann factors are stored this way. Since using arbitrary precision arithmetics slows the software down, the version using standard double-precision floating-point format it is also available.

Sorting flat structures

The Saffarian algorithm is specific by the fact that the size of P(i, j) for a state q U i,j or q P i,j is usually higher than for other DP schemes (Q, q root , ρ). The execution time of non-redundant sampling using Saffarian algorithm therefore greatly depends on the ordering of the flat structures within the list. Since a given p is chosen, for q U i,j or q P i,j from P(i, j) by generating the random number r and then summing up the probabilities p(λ | q, F S I , F S I next ) until the sum surpasses r the program will execute faster with lower number of structures. This can be achieved by ordering P(i, j), ∀, 1 ≤ i < j ≤ n by e -c p k B T by descending order. The effect of this sorting is shown on Figure 3.7; while in unsorted case we have to traverse seven different p, in second case, where they are sorted, it is only two. This becomes less effective for later stages of sampling since the low energy structures are picked first, however it considerably speeds up the initial stage of the algorithm. r r r r r r

Z χ Z χ Z χ Z χ Z χ Z χ Figure 3.8:
The effect of the sorting of list of available flat structures. During the random selection of the flat structure p from their set P(i, j), their sorting allows their faster selection. For the unsorted example (left), we need to traverse seven p before performing the selection, while we need to check only two in the case of sorted list (right).

The sorting of lists P(i, j), ∀, 1 ≤ i < j ≤ n was implemented into RNANR alongside the the non-redundant sampling algorithms.

Non-redundant sampling and the speed gain

The main motivation behind researching the non-redundant sampling is the faster generation of unique samples than for usual sampling approaches. We denote by U a set where any element appears at most once. A set of samples generated by non-redundant sampling has usually bigger coverage than the set of the same size with obtained using classic sampling.

Definition 3.3.1 (Coverage):

The coverage Cov(U) Cov(U) Cov(U) is is given by:

Cov(U) = S∈U Z S Z (3.14)
The coverage is therefore a proportion of partition function Z that is covered by the Boltzmann factors of U.

The bigger coverage implies the bigger probability that the next sample will be already in U unless sampled non-redundantly. In the case of non-redundant sampling, this is not a problem, since each sample is guaranteed to be unique. Consequently, non-redundant sampling saves time that would be otherwise spent by generation of duplicate secondary structures.

Definition 3.3.2 (Theoretical speed-up):

We call a theoretical speed-up a value E(|U|) gained for a sampling of a set of unique structures of size |U| by a non-redundant sampling:

E(|U|) = 1 + |U | i=0   1 - S∈U e -E S k B T Z   -1 |U| (3.15)
The theoretical speed-up is therefore an average number of times a single structure is redundantly sampled.

We computed theoretical speed-up on two sequences:

• 5S ribosomal RNA of Thermoplasma acidophilum (123 nt);

• Telomerase of RNA of Tetrahymena silvana (154 nt).

Results 3.3.1. The speed-up was computed using the Equation 3.15, with the se-quences being sampled and the necessary values being determined from these samples. The results indicate the exponential increase of the speedup as |U| increases (Figure 3.9). The length of the sequence has impact on the speed of increase and asymptotic limit -it is slower for longer sequences.

The important point of this example is that redundancy is an important factor when sampling the secondary structures and might noticeably slow down the sampling algorithms, justifying the research of non-redundant sampling methods. To compare the practical speed-up or gain of time, we reproduced the experiment from the paper introducing the locally optimal secondary structures sampling algorithm RNAlocmin based on scaling the sampling temperature [START_REF] Kucharik | Basin Hopping Graph: a computational framework to characterize RNA folding landscapes[END_REF] and using RNAsubopt from VIENNARNA as its basis. This experience compares

A ε = 5 % ε = 1 % ε = 0 . 1 %
RNAlocmin to another locally optimal secondary structure sampling algorithm -RNAlocopt [START_REF] William | Computing the partition function for kinetically trapped RNA secondary structures[END_REF]. All of the above programs use McCaskill algorithm with the DP scheme (L MC , q F 1,n , ρ MC ). The experiment consisted of tracing the number of unique secondary structures for certain time intervals. The sequence used is SV11, an artificially designed riboswitch generated by Qβ replicase [START_REF] Biebricher | Structural analysis of self-replicating RNA synthesized by Qβ replicase[END_REF]. We also compared the newer version of RNAlocmin. For the classical sampling algorithms -both versions RNAlocmin and RNAlocopt, we notice the exponential slow-down. This is related to the notion of theoretical speed-up -as the number of unique structures |U| increases, so does the probability of sampling an already existing secondary structure from U. The effect is specifically noticeable forRNAlocopt which does not handle the redundancy in any way, meaning for SV11 after certain time the new structures are basically not sampled. While RNAlocmin attenuates the effect of the redundancy to a certain degree by switching the sampling temperature, generating the new distribution of secondary structures where U has smaller proportion, the redundancy is still present.

On the other hand, RNANR does not sample the same structure twice and therefore the effect of the redundancy is not present, therefore the graph is linear and the sampling process by RNANR is, after short period of pre-computation notably of flat structures, much faster than by software employing classical approaches. Due to the sorting of flat structures it is expected that there should be slow attenuation of the speed once the considerable amount of secondary structures is sampled, however this does not pose a problem for usual sizes |U|. The practical result underline those obtained by analysis of theoretical speed-up -the redundancy is an important slow-down factor and its removal might speed up the sampling process considerably.

Landscape quality

While the previous analyses have shown that the non-redundant sampling is an important optimization factor for sampling DP based algorithms, it says nothing about the quality of the RNA landscapes that can be modeled from U. However, due to considerable lack of kinetic data regarding secondary structures of RNAs it is impossible to find a gold standard to which we could compare our results.

The second option is to pose a hypothesis on the properties of correct model of RNA folding landscape and perform the comparison between different software using a score based on these properties.

The good representation of a landscape should feature the most important key landmarks. This means not only potentially functional structures, but also the transitive states that represent the fastest folding pathways between said states. This implies that the best folding landscape should present the fastest folding kinetics, since slower kinetics implies missing some important transitive states. From this we can establish a measure of a quality of a folding landscape.

Definition 3.3.3 (Switching time):

Let P(S 1 , t) , resp. P(S 2 , t) be the proportions of sequences w formed a structure S 1 , resp. S 2 , and P(S 1 , 0) < P(S 2 , 0). The switching time t SW t SW t SW is the lowest time such that P(S 1 , t SW ) ≥ P(S 2 , t SW ).

Lower t SW implies faster kinetics and based on what precedes also the better landscape. Note that however if Z S 1 < Z S 2 , then it is possible that t SW → +∞ since the equilibrium implies higher Z S 1 . It is therefore important to choose correct states that will have their proportions studied.

The ideal candidates for the study would be riboswitches. Here the the first structure S 1 would be MFE (usually first functional structure) while the second would be the metastable structure -the second active state. This ensures t SW of finite time. The problem is the number of sequences of known riboswitches is limited.

For this reason, we created a set of 250 bistable sequences in a following manner:

• Create a sequence of length 100 nt with uniform distribution of nucleotides.

• Using RNANR, sample 1000 unique locally optimal secondary structures, taking parametrization into an account.

• Keep sequences where at least one secondary structure differing at least by 20 base pairs from MFE with them having less that 5kcal.mol -1 free energy difference. These structures are saved as well.

From these structures we generate a landscape using the different sampling software and the kinetic analysis described below using different software, which are then compared by their switching times t SW . Besides RNANR, RNAlocopt and RNAlocmin, we also compare RNASLOpt that detects stable locally optimal secondary structures [START_REF] Li | Finding stable local optimal RNA secondary structures[END_REF]. We also compare RNAsubopt from VIENNARNA library as a reference.

The pipeline for building an RNA landscape model from the set is based on the state-of-the art approaches detailed in Section 2.8. Since using barriers or BHG introduces too many transitive states and masks the results of the sampling by ameliorating them, we instead compute the activation energies and associated transition rates by using findpath and Metropolis rule respectively. In this case the landscapes are not represented by macrostates (see Section 2.9), but a network of locally optimal structures. The entire process of modeling an RNA folding landscape is as follows:

• By each software, sample a set of secondary structures U set of size 50, 75 and 100 samples.

• For each set remove duplicates and perform a gradient descent for each sample, generating U.

• From each U, create a graph G(V, E) where:

-V = U;
-E connects two structures S 1 and S 2 from U such that the number of different base pairs between them is less than K, where K is the lowest number for which G(V, E) is connected. The value of K is shared between all U made from U set .

• For each edge of G, compute the energy barrier using findpath, then associated transition rate by Metropolis rule (Section 2.8.2).

• Identify the two initial structures. Set the proportion of metastable state P(S meta , 0) = 1 and perform a kinetic simulation via numerical integration using treekin [START_REF] Michael | Efficient computation of RNA folding dynamics[END_REF]. Determine t SW .

Since the results are biased in the favor of RNANR that was used to establish the initial bistable sequences, we compensate this by giving a tolerance of difference of 10 base pairs for identifying the two initial states by the other software. The lowest t SW implies the best kinetics. ). The cases where not even RNANR retrieves both states implies that the state was not in first 50, 75 or 100 samples. However, the low retrieval percentage for other software also results from the redundancy, since after the removal of duplicates the number of samples in U was usually lower for them than for RNANR.

On the other hand, the switching time values t SW are consistently lowest for RNANR and with better proportion than successful identification of both metastable structures (Figure 3.11B). This implies that statistically RNANR returns better model of the RNA folding landscape under assumption that the faster kinetics implies better model. In the case that none of the software was successful at retrieval of both metastable structures, the sequence was marked as 'None', which is possible even for RNANR due to the sampling set used for bistable sequence generation having size of 1000. Any software failing to retrieve both structures was also disqualified. The results of the quality analysis by the different software. The successful identification shows RNANR was the most successful with retrieval of two metastable structures (A). This is expected due to RNANR being used to create the bistable sequences. However the lowest switching time is returned with better proportion than the retrieval of the structures (B). This hints at statistically better performance of RNANR when it comes to the analysis of RNA landscape. For (B), the 'None' category indicates the case where no software was able to retrieve both structures.

The concrete example of RNA folding landscape comparison is featured on Fig-ure 3.12 for he software RNANR and RNASLOpt, which returned the second best results after our software. For the given sequence w (Figure 3.12A) the bistable sequence presents an MFE structure (Figure 3.12D) and a metastable structure (Figure 3.12B). These structures are indicated on the landscapes generated from both software (Figure 3.12C). While in the case of RNANR the two structures are connected to two different halves of the landscape that are clearly weakly connected, in the case of the RNA landscape presented by RNASLOpt they are connected to the same half. This creates an auxiliary half in which w might get stuck due to limited transitions, making difficult to escape to state outside this part and slowing down the entire kinetics, which might be a reason the switching time between MFE and metastable structure is lower for RNANR. However, another reason for faster kinetics might be the higher connectivity of the folding landscape obtained from the samples generated by RNANR. Indeed, the higher connectivity makes the transition between the states easier, which may speed up the kinetics. In general, we noticed that the landscape models generated from the samples of RNANR exhibit higher edge to vertex ratio. To exclude the connectivity as a cause of faster kinetics, it would be good, in the future, to make analysis for the landscapes with the similar connectivity. RNANR both structures are connected to different parts while the landscape generated by RNASLOpt presents an auxiliary part. Likely for this reason the switching time is slightly lower for the landscape created by RNANR (E). However, the higher connectivity of the folding landscape generated for RNANR might also contribute to lower switching time.

The presented results suggest that RNANR and specifically non-redundant principle allow for better modeling of RNA folding landscapes, very probably due to an absence of redundancy and therefore presenting higher number of potentially interesting structures when constituting the landscape. For the future analyses, it may be interesting to study the landscapes with the similar connectivity for each software.

Non-redundancy for other DP schemes

While the first example employed Saffarian algorithm, using the introduced formalism we have shown that the non-redundant principle can be employed indeed for any DP scheme (Q, q root , ρ) as long as that scheme computes PfS. In theory it is enough to create an interface building a data structure d(w, Q, q root , ρ) for specific w that is adapted to said DP scheme and then connect the layer to the data structure of the DP scheme, an easy task due to number of interventions necessary for the DP scheme being rather limited.

To prove this point, we implemented the non-redundant sampling principle into

RNAsubopt from VIENNARNA library.

RNAsubopt

The RNAsubopt employs the DP scheme used in McCaskill algorithm

(Q MC , q F 1,n , ρ MC ).
The interface data-structure d(w, Q, q root , ρ) is built in a similar way than for the Saffarian algorithm, however the options for a given state q are not sorted like the flat structures and in plus the implementation also includes some options:

• Linking parent and children nodes:

-Hash table: the has function uses the values from q λ --→ c λ {q 1 , . . . , q u };

-Linked list: the same implementation as in RNANR (see Section 3.2.3);

• Numerical precision -this is the same as for RNANR:

-Double floating-point; -Arbitrary precision arithmetics: stores the partition functions of forbidden incomplete structures Z F S I in an arbitrary precision arithmetic format using MPFR; much slower but numerically stable when sampling a lot of structures or sequences with structures with important differences of free energy between them;

In this case the arbitrary precision-arithmetics was only included in the layer ensuring the non-redundant sampling, which was shown to be sufficient. That is convenient since it allows to limit the modifications of the DP scheme even more.

The performance profiling revealed that the memory allocation constitutes the major bottleneck of the non-redundant sampling program. For this reason the nodes of d(w, Q, q root , ρ) are not allocated alone, but in blocks that are preallocated, since this allows us to reduce the overhead associated to the single allocation. The size of the blocks is chosen sufficiently large to not allocate new blocks too frequently, but also sufficiently small to be easy to fit to available memory.

To establish a performance of the given implementation options, we extracted 365 sequences from the database RFAM [START_REF] Kalvari | Rfam 13.0: shifting to a genome-centric resource for noncoding RNA families[END_REF], specifically from the families RF00001, RF00005, RF00061, RF00174 RF01071 and RF01731 such that the %GC content of extracted sequences should be the most representative. For each sequence, we sampled 100, 200, 500, 1000, 2000, 5000, 10000, 20000, 35000 and 50000 samples by both the classical sampling and non-redundant sampling employing hash, linked lists and MPFR. We then analyzed the performance of the usual, redundant, and non-redundant versions of the sampling, both done by RNAsubopt from VIENNARNA.

The performance tests were executed on a laptop equipped with a CPU Intel R Core TM i7-5600U CPU with 2.60GHz × 4 on Linux Ubuntu 16.04 LTS and with 16 GB RAM.

Results

3.4.1.

Here we compare the obtained results for four representative selected sequences w of diverse lengths n:

• AE013598.1/1368227-1368452, from RF00174, length 226 nt;

• M84754.1/1-364, from RF00061, length 364 nt;

• Z50061.1/3-119 from RF00001, length 117 nt;

• AE015927.1/1705719-1705109 from RF01071, length 611 nt;

To shorten the notation, each w will be referred by the name of their family.

q q q q q q q q q q NR: 5. Execution time for sampling methods of sequence RF00174_0052 ( 226 nt ) q q q q q q q q q q NR: 3. Execution time for sampling methods of sequence RF00001_0012 ( 117 nt ) q q q q q q q q q q NR: 8. Execution time for sampling methods of sequence RF00061_0002 ( 364 nt ) q q q q q q q q q q NR: 0. At the first time, we consider the total number of unique structures sampled. We started by comparing of the version of the non-redundant sampling, implemented using hash table organization of the non-redundant layer constituted from the data-structure d, to the redundant sampling method. The results are shown on Figure 3.13. In general we notice the overhead of 50% independently of the length of the sequence, which is not satisfying. As stated, this is mostly due q q q q q q q q q q NR: 4. Execution time for sampling methods of sequence RF00174_0052 ( 226 nt ) q q q q q q q q q q NR: 2. Execution time for sampling methods of sequence RF00001_0012 ( 117 nt ) q q q q q q q q q q NR: 6. Execution time for sampling methods of sequence RF00061_0002 ( 364 nt ) q q q q q q q q q q NR: 0. were sampled for secondary structures by a redundant (blue) and non-redundant (light purple) method. The lines represent the linear approximation of the graphs. The time complexity overhead is in this case reduced to 10%-20%.

to the memory allocation, since nodes of the hash table are allocated separately.

For this reason we developed a second approach that implements the linked lists.

less objects than the hash table, reducing the overhead on object allocations. Second, the allocation was optimized by allocating a memory blocks instead of single objects. As a positive side-effect, the non-redundant sampling algorithm implementing the linked lists consumes less memory, allowing sample even longer sequences.

q q q q q q q q q q NR: 0. Execution time for sampling methods of sequence RF00174_0052 ( 226 nt ) q q q q q q q q q q NR: 5. As a side note, we also provide a reduced example comparing the efficiency of the non-redundant sampling implemented via linked lists where the non-redundant sampling layer stores the values using the MPFR library. The results for the two sequences from the families RF00174 and RF00061 can be consulted on Figure 3.15. Using the arbitrary precision arithmetic is accompanied by the considerable augmentation of the time complexity, in the case of the sequence from the family RF00174 by rising from 10% to 190%. Therefore, while MPFR may be necessary to use when sampling a more delicate sequences, such as the artificially designed riboswitch SV11 generated by Qβ replicase [START_REF] Biebricher | Structural analysis of self-replicating RNA synthesized by Qβ replicase[END_REF] which present considerable free energy differences between the secondary structures of interest, it is best to avoid it.

We now concentrate on the performance of both sampling methods in regards to the number of unique secondary structures sampled. For the classical, redundant q q q q q q q q q q q q q q q q q q q q 0 10000 20000 30000 40000 50000 0.0 0.5 1.0 1.5 2.0 Execution time (s) Number of samples Sampling Types q q q q R sampling NR sampling

Sampling speed of sequence RF00174_0052 ( 226 nt ) q q q q q q q q q q q q q q q q q q q q 0 10000 20000 30000 40000 50000 0.0 0.5 1.0

Execution time (s)

Number of samples Sampling Types q q q q R sampling NR sampling Sampling speed of sequence RF00001_0012 ( 117 nt ) q q q q q q q q q q q q q q q q q q q q 0 10000 20000 30000 40000 50000

0 1 2 3
Execution time (s) Number of samples Sampling Types q q q q R sampling NR sampling

Sampling speed of sequence RF00061_0002 ( 364 nt ) q q q q q q q q q q q q q q q q q q q q 0 10000 20000 For the shorter sequences A and C, the non-redundant sampling performs better due to high number of the duplicates generated by the usual sampling approach. For the sequences B and D, it performs worse due to the complexity overhead and little probability to generate duplicates for such long sequences at the beginning.

sampling, we remove the duplicates, and then we compare the number of the secondary structures sampled at the same time. The results are show on Figure 3.16.

For the shorter sequences from families RF00174 and RF00001 the non-redundant sampling gives higher number of unique structures for the same execution time than the usual sampling approach. This is because for the shorter sequences it is probable to generate a duplicate quite early on, therefore the usual sampling methods tend to generate lot of them. On the other hand, this is unlikely to happen for long sequences early on during the sampling, since there is much more options to generate a secondary structure, and consequently the probability of generating an already sampled structure is low. In this case, it is the complexity overhead that slows the non-redundant sampling estimator when compared to its classic counterpart, consequently performing worse.

We have demonstrated in Section 3.3.3 that the speed gain becomes more important the more structures we sample. For the longer sequences the slow down from the redundancy eventually causes the performance of the usual sampling method to fall enough so that of the non-redundant sampling overcomes it, but it may take a considerable number of secondary structures to sample. The extensive, long performance test with high number of samples is one of the analyses that may be interesting to perform in the future.

Conclusion -Non-redundant sampling

In this section, we presented the concept of non-redundant sampling. The formalism that we introduced in the previous chapter allowed us to demonstrate how its principle can be implemented into all DP schemes that compute the partition function. The non-redundant sampling necessitates a creation of a layer with specific data-structure to function, notably it presents a way to store the values associated to the sets of structures that still can be generated. However, we have demonstrated that the layer has limited number of interactions with the DP scheme itself, making it easy to implement. Moreover, the numerical instability issues can be solved by implementing the arbitrary precision arithmetics uniquely to the data-structure related to the non-redundant sampling, limiting the modification of the DP scheme layer even more.

The non-redundant sampling was tested for the Saffarian algorithm [START_REF] Saffarian | RNA locally optimal secondary structures[END_REF], implemented in RNANR jointly with H. Touzet and the McCaskill algorithm used by RNAsubopt from VIENNARNA library [START_REF] Lorenz | Vi-ennaRNA Package 2.0[END_REF]. Both software show the some advantages over the competing software. Notably, RNANR samples unique secondary structures more efficiently than its competitors and it seems to return better predictions of RNA folding landscapes. However, it may be interesting, in the future, to study and compare the folding landscapes with similar connectivity. The non-redundant sampling implemented in RNAsubopt shows a better perfor-mance regarding sampling unique secondary structures for shorter sequences.

The optimized implementation is an important part of the efficiency of the nonredundant sampling. While the conception of the non-redundant sampling layer itself is easy enough, to function efficiently it must be as optimized as possible. This section has shown the impact of the correct implementation of the layer, where the replacement of the hash tables by linked lists reduced the complexity overhead by 30%. It is also important to note that using the arbitrary precision arithmetics considerably slows the algorithm down, so it is best to use it only when really needed.

One caveat is that for sampling a limited set of secondary structures for long sequences by non-redundant sampling is inefficient. This is mainly due to complexity overhead, which may be limited but not completely removed. However, this effect is expected to disappear once a sampling set of sufficient size was sampled, since the effect of the redundancy amplifies with it.

Chapter 4

Estimating from non-redundant sampling sets

The non-redundant sampling samples the unique secondary structures. While this means faster coverage of the associated folding landscape, it is also accompanied by the introduction of a dependency between the samples as well as lacking the information about their frequency. This proves problematic when one wants to estimate a certain quantity of such a sample.

The absence of the frequency means that the usual estimator cannot be employed. However, we have researched, the estimator specific to non-redundant samples. This estimator proves better convergence than the usual approach employed for classic sampling sets, and exhibits smaller variance in most cases.

In this section, after introducing the necessary notions and the classic way of estimating a certain quantity for redundant sampling sets, we explain how to obtain the equivalent for the non-redundant sampling sets. This will be followed by experiments using the implementation of non-redundant sampling for VIEN-NARNA library on the sequences from the RFAM database [START_REF] Kalvari | Rfam 13.0: shifting to a genome-centric resource for noncoding RNA families[END_REF] cases to estimate a number of representative quantities and the the analysis of the influence of different properties of investigates sequences on the quality of the estimation.

The work described in this whole section was done in collaboration with Christelle Rovetta, post-doctoral fellow within the AMIBio team.

The non-redundant estimator

The estimators in general have for objective to give an estimation of a certain quantity, represented by a specific function.

F(S) =    1 if (i, j) ∈ S 0 otherwise (4.1)
We can therefore specify the central problematics of this section. • Input: A sample set U drawn from non-redundant distributions with probability of S ∈ U:

p(S | F) = Z S -1 S∈F × Z S Z -Z F . (4.2)
• Output: An unbiased estimation F(U) of F for U equal to the expectation

E(F(U)).

There are many usual approaches to estimate the probability of the classic, redundant samples, however, for the non-redundant generation the situation is a bit more complicated.

Empirical mean

Let X be the classical, redundant sampling set of ordered samples S 1 , . . . , S |X | and

X i = (S 1 , . . . , S i ) for 1 ≤ i ≤ |X |.
The empirical mean is given by

F(X ) = 1 |X | S∈X F(S).
This estimator uses the frequency of the samples as a way to approximate its probability. This means that the estimator is unbiased, because

E (F(X )) = 1 |X | S∈X E (F(S)) = E (F(S n )) .
As a reminder, S n is the set of all secondary structures for w of length n, constituting a variable with |S n | states.

In the case of non-redundant sampling set U the above no longer applies, as the absence of frequencies means the probability cannot be expressed by it. We need a way to express the probability for U before the non-redundant sampling set specific estimator can be established.

Example 4.1.2. Consider 100 identical RNA secondary structures with a single base pair and 50 identical molecules with two base pairs. The average number N BP of the base pairs in such set is:

E(N BP ) = 1 150 100 i=1 2 + 50 i=1 1 = 4 3 (4.3)
But if we try to use such approach on exhaustive non-redundant sampling set, of size 2, we get:

N BP = 1 2 × (1 + 2) = 3 2 (4.4)
which is not equal to the expected value E(N BP ) = 4 3 .

Removing the dependency from the estimator

We remind that the probability p(S | U) of non-redundantly choosing a sample S as:

p(S | U) = Z S -1 S∈U × Z S Z -Z U . (4.5)
In this section, we suppose U ⊆ S n is ordered, S 1 , • • • , S |U | ∈ U being generated in this order. We can define

U i = ( S 1 , • • • , S i ) with 1 ≤ i ≤ |U| as a subset of U.
First we establish a k th moment of U.

Definition 4.1.2 (k th moment of set): he k th moment of a set U is the value:

µ (k) U (F) = S∈U p( S | U) × F( S) k
With this, we can define the non-redundant estimator. . Replacing the probability by its expression from Equation 4.5 gives:

F(U) = 1 |U|   |U | i=0 F( S i ) 1 -µ (0) U i-1 (F) + µ (1)
U i-1 (F)   .
Since E N BP ( S 1 ) = E N BP ( S 2 ) = E(N BP ) = 4 3 , this becomes E N BP (U) = 1 2 4 3 + 4 3 - 4 3 × p( S 1 | U)) + p( S 1 | U) × 4 3 E N BP (U) = 1 2 8 3 E N BP (U) =
E U (F( S)) = S∈Sn Z S -1 S∈U × Z S Z -Z U × F( S) E U (F( S)) = S∈Sn/U Z S Z -Z U × F( S) since if S ∈ U then p( S | U) = 0. It follows E U (F( S)) = S∈Sn/U p(S) 1 -µ (0) U (F) × F( S) E U (F( S)) = 1 1 -µ (0) U (F)   S∈Sn/U p(S) × F( S)   E U (F( S)) = 1 1 -µ (0) U (F)   S∈Sn p( S) × F( S) - S∈U p( S) × F( S)   .
since U ∈ S n . The term S∈Sn p( S) × F( S) = E(F(S n )) since it is estimated on S n . This gives us

E U (F( S)) = 1 1 -µ (0) U (F) E(F( S)) -µ (1)
U (F)

E(F( S)) = E U (F( S)) × (1 -µ (0) U (F)) + µ (1)
U (F).

The right-hand side of the expression can be used to establish an estimator F(U)

F(U) = 1 |U|   |U | i=0 F( S i ) 1 -µ (0) U i-1 (F) + µ (1)
U i-1 (F)   . (4.6)
It remains to prove that the estimator F(U) is unbiased, meaning we have to show that

E F(U) = E (F(S n )) .

123

This can be done by induction on the size of U. Base case (|U| = 1): In this case U 0 = ∅ and µ (0)

U 0 = µ (1)
U 0 = 0. From that we have:

E F(U) = E F( S 1 ) = E (F(S n ))
Induction (|U| = m + 1): Suppose that the property is valid for |U| = m and we have to verify it for |U| = m + 1. We have

E F(U) = E (F(S n )) E F(U) = E 1 m + 1 m+1 i=0 F( S i ) 1 -µ (0) U i-1 (F) + µ (1) U i-1 (F) E F(U) = E m m + 1 F(U m ) + 1 m + 1 E F( S m+1 ) 1 -µ (0) Um (F) + µ (1)
Um (F)

E F(U) = m m + 1 E (F(S n )) + 1 m + 1 E F( S m+1 ) 1 -µ (0) Um (F) + µ (1)
Um (F) .

The only thing left is to prove that E F( S m+1 ) 1 -µ (0)

Um (F) + µ (1)
Um (F) = E (F(S n )). We pose B m+1 = F( S m+1 ) 1 -µ (0) Um (F) + µ (1)
Um (F).

We define p(U m ) as a probability of having a set of unique structures U m . We then develop E(B m ):

E(B m ) = Um||Um|=m p(U m )E Um F( S m+1 )(1 -µ (0) Um (F)) + µ (1)
Um (F)

E(B m ) = Um||Um|=m p(U m )E Um   µ (1)
Um (F) + S∈Sn/Um p( S | U m )F( S)(1 -µ (0) Um (F))   E(B m ) = Um||Um|=m p(U m )   µ (1)
Um (F) + S∈Sn/Um p( S) 1 -µ (0) Um (F) F( S)(1 -µ (0) Um (F))   E(B m ) = Um||Um|=m p(U m )   µ (1)
Um (F) +

Um (F) + E(F(S n )) -µ (1)
Um (F)

E(B m ) = Um||Um|=m p(U m )E(F(S n )) E(B m ) = E(F(S n ))
from which it follows directly that

E F(U) = E(F(S n )).
The proposed estimator however also presents one important disadvantage of it being inconsistent. When the number of the employed sampling set raises indefinitely, it does not converge towards the true value. While the estimator is unbiased, meaning its expected value F(U) gives the true value of F(S n ), it does not converge towards it for the entire sampling run. This stems from the fact that the number of non-redundant samples that can be generated is limited. This issue is addressed at the end of this Section.

When it comes to variance, the numerical analysis does not allow us to conclude decisively on its value when comparing F(U) and F(X ). In practice it also depends on the complexity overhead of the implementation, as the same-size nonredundant set U cannot be produced in the same time as the classic set X due to the additional modifications to be made by the non-redundant sampling layer. The experiments presented in the next section compare the performance of both estimators.

Applications of non-redundant estimator

Here we present a number of applications of non-redundant estimator to specific cases and the comparison of the performance of the classic estimator F(X ) applied on sampling set X and that of the classic, redundant estimator F(U). For every experiment, unless noted otherwise, we used the dataset extracted from RFAM database [START_REF] Kalvari | Rfam 13.0: shifting to a genome-centric resource for noncoding RNA families[END_REF], specifically 365 sequences of RFAM families RF 00001, RF 00005, RF 00061, RF 00174, RF 01071 and RF 01731. The sequences were chosen in a manner so that their GC content is as homogeneous as possible. Both the redundant, from this point on referred as classical, and non-redundant sampling is performed by RNAsubopt of VIENNARNA library, the latter of which was implemented by us as discussed in Section 3.4.1.

Estimating base pair probabilities Definition (Dotplot matrix):

A dotplot matrix is a matrix D associated to a sequence w storing the values:

D i,j = S∈Sn    p(S) if (i, j) ∈ S 0 otherwise
In other words, D i,j is the average probability of the given structure having a base pair (i, j). This matrix can be computed in an exact manner by RNAfold from VIENNARNA library. Therefore, we compute the reference value D using RNAfold and the estimations by using the classic mean average estimator D(X ) and the non-redundant sampling specific estimator U introduced by us. For both U and X , we sampled 10 3 structures. We then computed their respective errors: ; for X the unique set was generated from them), length n of w and the GC-content. the difference between e R and e NR , and the color the estimator that performs better (blue if e R > e NR , orange otherwise). There does not seem to be an influence of the GC-content, though D(U) seems to perform better for the sequences where GC-content is biased (Figure 4.1A). This is a little surprising, considering that higher GC-content would imply stronger structures due to higher number of G-C bonds, bigger free energy differences and consequently more redundancy.

On the other hand, there is an influence of the n as well as Cov (Figure 4.1B and C resp.). The estimator D(U) performs better than D(X ) for shorter sequences and higher Cov. This is logical since the complexity overhead is not as noticeable for shorter sequences. Similarly, the effect is more remarkable the bigger percentage of the search space is covered.

We also include the Figure 4.2 that indicates the total number of cases where U performs better. While the difference between the errors is in all cases close to 0, the non-redundant estimator D(U) performs better in around 300 cases out of 365. Coherently the majority of cases where D(U) performs better is those with high coverage Cov, the intervals of which are marked by the color code. Due to inconsistency, the number of sampled structures must stay relatively small. The example of dotplot matrix D indicate that in the majority of the cases the estimator specific to non-redundant samples performs better. On the other hand, the difference of performances is smaller for longer sequences, which, despite bing expected and coherent with the observations from Section 3.4.1 where the performance of the non-redundant sampling was worse for long sequences and small |U|, is still a little disappointing since the main reason behind the establishment of the non-redundant sampling is the analysis of the longer sequences.

For two specific sequences w, M30199.1/68 -167 from the family RF00001 and CP000679.1/1996671 -1997302 from RF01071, which will be referred by their specific families, we also compared the overall performance of both estimators D(X ) and D(U). with the respect to the number of performed samples and the execution time. The maximum number of sampled structures for both w and sampling sets U and X was 10 4 . For the time performance, we sampled the sampling sets U and X for the approximately same execution time. The sampling was executed on a laptop equipped with a CPU Intel R Core TM i7-5600U CPU with 2.60GHz × 4 on Linux Ubuntu 16.04 LTS and with 16 GB RAM. We then analyzed the evolution of e R and e NR in function of the number of samples and execution time. Also, notice that in the case of the sequence from the family RF01071, the error starts to raise after some time, resp. number of samples (Figure 4.3 left bottom). This might point to an inconsistency of our estimator. While our estimator is unbiased, meaning its estimated value is equal to the expected value of given feature function, here D, it does not converge towards the predicted value for an infinitely large sample. This problem is addressed later.

Graph distance distribution

Here we suppose the RNA secondary structure S as a graph G S = (V, E) where:

• V are nucleotides of the sequence;

• E connect two nucleotides connected by the backbone or a base pair. Here we first test a distribution of d i,j (S) ∀S ∈ S n , denoted D i,j (S) generated by the sampling methods themselves. Unlike for the dotplot matrix application, there is not an implemented and efficient way to compute the minimum graph distance distribution for secondary structures in an exact manner. However, we may exploit the consequences of the Chebyshev's inequality. Theorem 4.2.1 (Chebyshev's inequality): The Chebyshev's inequality is, for a random variable X with finite expected value µ and a finite non-zero standard deviation σ and a real positive number k, given by:

p((X -µ) ≥ kσ) ≤ 1 k 2
In the other words, only a part of the values X will be more than a certain dis-tance of the mean µ. The consequence of this inequality is the fact that for sufficient number of the samples a significant proportion of them will be distributed around the value µ. We can therefore generate the reference structure by sampling redundantly sufficiently extensive number of the secondary structures and computing D i,j (X ).

In our case the population is all secondary structures S ∈ S n . By exploiting the Chebyshev's inequality we compute the reference distribution D i,j (S n ), then compare it to the results of both classical sampling d i,j (X ) and the non-redundant one D i,j (U). Due to the extensive computations necessary to generate a reference, this experiment is performed only for the 100 nt long sequence M30199. The non-redundant sampling distribution D 16,82 (U) is closer to the reference, proving the non-redundant sampling to be more efficient in this case. Also note that U shows more distance classes that X and is closer to the expected value -this is not surprising since non-redundant sampling finds more unique structure and consequently has better chance to find more of them.

To test the estimators of the both sampling approaches, we use a feature function of the weighted graph distance of all possible couples of nucleotides. Again, M is computed from 10 6 samples produced by classical sampling by exploiting the consequences of Chebyshev's inequality and is taken as a reference. Both e R ande NR were computed for up to 10 4 samples.

For the previous experiment, we also determined two values B 1 and B 2 :

• B 1 is the number of samples necessary to have e R > e NR for the first time;

• B 2 is the number of samples necessary to always have e R > e NR for any subsequent samples.

The above experiment was repeated 50 times, each producing each iteration that gave us specific X and U, the values B 1 and B 2 were calculated.

We then plotted their distribution. For the analysis of the values of B 1 and alwaysbest, we notice that the non-redundant estimator shows, for the sequence M30199.1/68-167, the better value for the first time after sampling about 2500 samples at average. The value becomes consistently better after sampling close to 3200 structures. Interestingly, the non-redundant estimator shows at average better values that the classical estimator despite being inconsistent, though this can be explained by not sampling sufficiently big number of the secondary structures. The experiments on minimum graph distance and graph distance matrices notably confirm that the estimator F for non-redundant sampling sets U has smaller variance that F for X , as well as giving the estimation that for the reasonable number of secondary structures are closer to the true value.

Shape probabilities

In this experiment, by RNA shape we denote a simplification of dot-bracket format for RNA secondary structures [START_REF] Giegerich | Abstract shapes of rna[END_REF]. There are many ways how one can simplify a secondary structure, so in our case we will define it in a following way. Definition 4.2.4 (Shape): The shape abstraction mapping π [START_REF] Giegerich | Abstract shapes of rna[END_REF] is a transformation f SHA (S) of a dot-bracket format of a secondary structure S where:

• Dot '.' is removed.

• Helix is replaced by a single pair of square brackets; The SHAPE is the last example on which the estimators were tested. We define a feature function F SHA such that

F SHA (S n ) = S∈Sn    1 if f SHA (S) = f SHA (S MFE ) 0 otherwise
where S MFE is the MFE structure. In other words, F SHA (S n ) is the probability that the shape of S ∈ S n is identical to MFE.

The value of F SHA (S n ) is estimated by the non-redundant estimator F(U) for a set of unique samples U and by its classical equivalent F(X ) for redundant sampling set X . This value was evaluated for the sequences w and respective families:

• M30199.1/68 -167 from RF00001;

• AY344021.1/1 -348 from RF00061;

• CP000283.1/2593935 -2594143 from RF00174;

• CP000679.1/1996671 -1997302 from RF01071.

To simplify the notation, we will refer to these sequences by the name of their respective families.

Like for the previous application, the computation of the exact value of F SHA (S n ) is difficult, therefore we resort to estimating the value from sufficiently large sampling set of 10 6 . For each w, we sampled |X | = |U| = 10 4 secondary structures.

Resolving the convergence problems

The inconsistency related to the non-redundant sampling estimator F(U) pertains to the fact that non-redundant sampling has limited number of samples that can be produced. We do not dispose with the existing values of F( S j ) for any j > l where l is the maximum number of the unique structures, resulting in the inconsistency. This short section presents a way to make the theoretical value of the estimator converge in such cases.

Let us define ∀j > l, F( S j ) = 0, and m > l. In that case, the non-redundant estimator

F(U) = 1 m m i=0 F( S i ) 1 -µ (0) U i-1 (F) + µ (1)
U i-1 (F) can be divided in two separate terms

F( U) = 1 m l i=0 F( S i ) 1 -µ (0) U i-1 (F) + µ (1)
U i-1 (F) + m i=l+1 F( S i ) 1 -µ (0) U i-1 (F) + µ (1)
U i-1 (F) .

For the second term, we have

m i=l+1 F( S i ) 1 -µ (0) U i-1 (F) + µ (1)
U i-1 (F) = (m -l) × E(F(S n ))
since the first term is nullified by F( S i ) = 0 and for the second we have µ

U (F) = S∈U p( S | U) × F( S) (1)
which in this case completely covers the space, and therefore is equal to E(F(S n )).

Consequently, we have

F(U) = 1 m F(U l ) + (m -l) × E(F(S n )) F(U) = 1 m × F(U l ) + m -l m × E(F(S n )). (4.7)
When we test the solution for the convergence when m → +∞, we find

lim m→+∞ F(U) = E(F(S n )) (4.8)
since the first term of Equation 4.7 becomes negligible and the second becomes close to 1 for values of m sufficiently large. Therefore, substituting F( S j ) by 0 for j > l makes the non-redundant sampling estimator indeed converge towards the expected value and therefore the estimator becomes consistent.

Conclusion -Non-redundant estimator

In this section we discussed the problem of the estimating a specific quantities or properties, represented by a quantifiable feature function, of non-redundant sampling set. The lack of the frequencies forbids us from simply using the weighted average as non-redundant estimator. On the other hand, using the conditional probability leads us to an unbiased estimator that allows us to estimate these quantities from non-redundant sampling.

The advantage of this estimator seems to be its smaller variance, as shown by the practical applications on different quantities such as base pair probability or minimum graph distance. In the former case, our estimator performed worse in about 65 cases, having an important edge over the classical sampling approach and the quantity estimated by the weighted average. The better performance was noted specifically for shorter sequences and bigger coverage, however the non-redundant estimator still seems to perform better for long sequences than the classic counterpart, even if with a smaller difference. It is also important to notice that our estimator can perform worse in some specific cases, such as small sampling sets for long sequences where the redundancy does not really matter.

The fact that the non-redundant estimator performs better than the classical version not only with the respect to the number of samples, but also the execution time demonstrates that the non-redundant sampling can be implemented efficiently into DP sampling schemes, with little complexity overhead. However, it may need some considerable optimization, like was the use of the linked lists for the implementation used in RNAsubopt in VIENNARNA.

The one problem related to the non-redundant sampling set estimator is the fact that the estimator is inconsistent, or it does not converge towards its expected value even when it is unbiased once the number of samples becomes sufficiently big. This is because the non-redundant sampling set is of limited size. For theoretical computations, we solved this problem by substituting the value of feature function of all theoretical samples exceeding the number of possible unique structures by 0, which makes the estimator convergent.

Chapter 5

Extending towards pseudoknots

Until now, we assumed, as stated in Section 2.1, that the secondary structures contain only the base pairs that cannot cross each another. This condition is necessary to make the decomposition into loops, the basis of all the DP schemes detailed in this work possible since we can assume the independence between what happens inside a given base pair and what happens outside of it. This hypothesis also allows the computation of the free energies of the secondary structures by the Turner Energy Model.

However, such simplification removes an entire class of substructures from the equation -the pseudoknots. We already stated in Definition 2.4.2 that the pseudoknot consists of crossing base pairs. There are many different types of pseudoknots from a simple crossing of two helices to complicated structures [START_REF] Christian M Reidys | Topology and prediction of RNA pseudoknots[END_REF]. While complex pseudoknots are unlikely to appear in biologically important structures, many simpler pseudoknots have important role in biology [START_REF] Theimer | Structure of the human telomerase RNA pseudoknot reveals conserved tertiary interactions essential for function[END_REF][START_REF] Jacobs | Identification of functional, endogenous programmed-1 ribosomal frameshift signals in the genome of Saccharomyces cerevisiae[END_REF].

The reason for such restriction is simple -reduce the time complexity. It was demonstrated that accounting for a pseudoknot in its most general definition is an NP-complete problem [START_REF] Akutsu | Dynamic programming algorithms for RNA secondary structure prediction with pseudoknots[END_REF][START_REF] Rune | Pseudoknots in RNA secondary structures[END_REF]. Even the most efficient algorithms for the simplest classes of pseudoknots have O(n 4 ) time complexity compared to O(n 3 ) of the usual DP schemes. However, due to their biological significance, it may be interesting to analyze the simpler pseudoknots and to propose an efficient sampling method on which a non-redundant sampling principle can be eventually extended.

The work in this section is preliminary and still in progress. For this reason, this section is more the presentation of the concept and the main ideas than the re-sults and the completed software. We begin by the introduction of the algorithm behind the software pknotsRG [START_REF] Reeder | pknotsRG: RNA pseudoknot folding including near-optimal structures and sliding windows[END_REF] that was also included in pKiss [START_REF] Janssen | The RNA shapes studio[END_REF]. Both software predict the structures with the most optimal pseudoknots for a given interval. The one returns pseudoknots consisting of two perfect helices crossing, while the second returns also kissing hairpins, where one helix crosses once with two others. We will then discuss how to generalize this algorithm while keeping the complexity reasonable, and how to extend it on suboptimal structures. Finally, we will discuss how to extend the non-redundant sampling principle on proposed algorithm and on the future work that needs to be done in this direction.

State of the Art -pknotsRG algorithm

Here we explain the algorithm used in pknotsRG [START_REF] Reeder | pknotsRG: RNA pseudoknot folding including near-optimal structures and sliding windows[END_REF]. This program predicts the so-called canonical simple recursive pseudoknots.

Definition 5.1.1 (Canonical simple recursive pseudoknot): A canonical simple recursive pseudoknot consists of exactly two helices H 1 (i, j, lh) and H 2 (i , j , lh ) with i < i , of maximum possible length lh and lh and such that i + lh < i < i + lh < j -lh < j < j -lh.

The helices, as defined in Definition 2.4.2, must not contain an unpaired base pair. Such helices are also called perfect helices. The simple recursive pseudoknot is then a crossing of two such helices. The maximality of lh, lh implies that the helix is extended until it is not possible to create a base pair anymore. If the helices would overlap, H 1 (i, j, lh) with lower i takes the priority and is extended to its maximum and H 2 (i , j , lh) takes the remaining base pairs. This also completely defines the internal pseudoknot segments w[i + lh, i -1], w[i + lh, j -lh] and w[j + 1, j -lh] which we respectively denote I 1 , I 2 and I 3 . The example of a canonical simple recursive pseudoknot is given on Figure 5.1.

The algorithm to search for said pseudoknots is designed in a following manner:

• For each substring w[i, j] with i < j, pre-compute and store the maximum length a of a perfect helix H(i, j, a) in a look-up table.

i i i i + lh -1 i + lh -1 i + lh -1 j j j j -lh + 1 j -lh + 1 j -lh + 1 i i i i + lh -1 i + lh -1 i + lh -1 j j j j -lh + 1 j -lh + 1 j -lh + 1 I 1 I 1 I 1 I 2 I 2 I 2 I 3 I 3 I 3 H 1 (i, j, lh) H 2 (i , j , lh ) a) b) i i i i + lh -1 i + lh -1 i + lh -1 j j j j - l h + 1 j - l h + 1 j - l h + 1 i i i i + l h - 1 i + l h - 1 i + l h - 1 j j j j -lh + 1 j -lh + 1 j -lh + 1 H 1 (i, j, lh)
H 2 (i , j , lh ) • For each quadruplet (i, j, i , j ), if H 1 (i, j, a) and H 2 (i , j , a ) do not overlap, meaning i + a -1 < j -a + 1, create the pseudoknot directly. If they do, set a = j -i -a + 1 and create the pseudoknot. Compute its energy and for each w[i, j], keep only the pseudoknot with the lowest one.

I 1 I 1 I 1 I 2 I 2 I 2 I 3 I 3 I 3
The complexity of this algorithm is O(n 4 ) in time and O(n 2 ) in space. The limiting factor is to analyze all quadruplets (i, j, i , j ), giving the complexity of O(n 4 ) since once all indexes are known the the pseudoknot is assembled from the values in the look-up table which are computed in O(n 2 ) complexity. In the base version, only the most optimal pseudoknot is stored for given substring w[i, j], and the pre-computation step stores the maximum length of helix starting with (i, j), giving the memory complexity of O(n 2 ).

The free energy of a pseudoknot is computed as the sum of the energies of H 1 (i, j, a) and H 2 (i , j , a ) and adding a penalty for creating a pseudoknot. The pseudoknot penalty was not determined empirically and therefore may be inaccurate.

This algorithm is employed by both pknotsRG and pKiss. Additionally, pKiss allows to search for canonical simple recursive kissing hairpins by four different methods, which are not discussed here. The webservice pKiss can be consulted at https://bibiserv2.cebitec.uni-bielefeld.de/pkiss.

There are many other algorithms that compute the pseudoknots of varying complexity and specificity. However, many of them are too stringent or too general for our application, while the algorithm discussed here has low complexity and needs only to be generalized to encompass a wider class of pseudoknotted structures.

Generalizing the pknotsRG algorithm

The algorithm employed in pknotsRG as well as pKiss is useful for discovering the simplest of pseudoknots in an efficient manner. However the number of returned pseudoknots is restricted and due to choosing only the best pseudoknot for given w[i, j], it does not propose suboptimal solutions for said structure. Therefore, we want to propose an extension of the previous algorithm.

Problem 5.2.1:

• Input: A sequence w;

• Output: Set of secondary structures S that contain pseudoknots produced by a combing two, not necessarily perfect, helices in a Boltzmann distribution.

As said, there are many existing algorithms that can generate pseudoknotted structures; the algorithm of pknotsRG is only one example that was chosen as a starting point due to its complexity and implementation advantages. Other examples would be Akutsu algorithm [START_REF] Akutsu | Dynamic programming algorithms for RNA secondary structure prediction with pseudoknots[END_REF], which likewise treats only pseudoknots consisting from two helices, or that established by Elena Rivas and Sean R. Eddy [START_REF] Rivas | A dynamic programming algorithm for RNA structure prediction including pseudoknots[END_REF], capable to predict even more complicated pseudoknots but at higher time and memory costs. Our algorithm, while not as fast as that used in pknotsRG, is able to predict multiple pseudknotted secondary structures, including those with suboptimal pseudoknots, that pknotsRG is unable to find, justifying its research.

We first introduce the concept of imperfect helices. We then continue how to build a list of helices, which the may be used to assemble the pseudoknots.

Imperfect helices

Our first objective is to establish a table of all helices that can be used to specifically build the pseudoknots we are interested in. In our case these helices may be imperfect.

Definition 5.2.1 (Imperfect helix): An imperfect helix h(i, j, k, l) is a set of base pairs constituting the helix and being delimited by (i, j) and (k, l). w[i, k] and w[l, j] may contain unpaired bases.

n∆ of them due to the position of l being determined by ∆ and that of k. This effectively reduces the time complexity of searching for helices to O(n 3 ∆). The memory complexity is likewise reduced.

The introduction of the H-type pseudoknots K H ( H 1 , H 2 ) forces us to redefine the space of all secondary structures S n , since the condition of no-crossing base pairs is no longer valid. From this point on, we suppose S n contains all secondary structures S ∈ S n that may contain elements K H ( H 1 , H 2 ) with ∆ in effect. All other non-conflicting restrictions still apply as defined in Section 2.1.

We are now interested to enumerate all imperfect helices starting with a base pair (i, j). For this reason, we introduces a notation of a set of imperfect branches

H ∆ i,j H ∆ i,j
H ∆ i,j a set of all imperfect branches H(i, j, k, l) with i ≤ k < l ≤ j containing only helices h(i, j, k, l) that validate ∆.

To keep the complexity of the algorithm building the helices reasonable, we propose the following dynamic programming scheme. Each base pair (i, j) can be either itself alone a helix H(i, j, i, j) or it may extend some H(i , j , k, l) such that i < i < j < i . We also need to know the partition function Z H of each helix. To simplify a notations, let Z H = Z H(i ,j ,k,l) . Therefore:

H ∆ i,j =              (i, j, 1) |j-j -(i -i)|<∆ |j-l-(k-i)|<∆ j-j <∆ i -i<∆ {(k, l, Z H × e E ILG (i,j,i j ) k B T ) | (k, l, Z H ) ∈ H ∆ i ,j } . (5.1)
We can then enumerate all possible helices in a recursive manner.

The conditions j -j < ∆ and i -i < ∆ come from the restriction of the length of unpaired nucleotides in h(i, j, k, l) ∈ H(i, j, k, l). This way we can enumerate all possible imperfect helices starting with (i, j), including suboptimal helices.

Originally, H ∆ i,j was intended to be stored in form of posets, see Section 5.2.3. Since the performance of such a solution was worse than expected, we use a simple table of n × (∆ + 1) dimensions that stores the sum of Z H of all helices H(i, j, k, j + i -k -δ), where δ = j -l -(k -i) and therefore |δ| < ∆. This way we can look directly whether any imperfect helix h(i, j, k, l) exists or not. The imperfect exact helix can be obtained by performing the stochastic backtrack on the value Z from the position (i, j, k, j + i -k -δ) in O(n × ∆) time if i, j are known. While this makes the values of H ∆ i,j not necessary, the relation stays very useful when pre-computing all helices for w[i, j].

Assembling the pseudoknots

The generation of all possible pseudoknots K H ( H 1 , H 2 ) for a given substring w[i, j] is rather simple once we pre-compute the table of imperfect helices. We need to simply adapt the previous DP scheme. Definition 5.2.4 (Pseudoknot space): We call a pseudoknot space K (i, j), or an H-type pseudoknot space, a space of all pseudoknots spanning over substring w[i, j] defined as

K (i, j) = k<i <k <l j <l { H 1 (i, j , k, l) ∪ H 2 (i , j, k , l )}.
The internal segments I 1 , I 2 and I 3 are not considered.

The space K (i, j) can be generated as follows:

• Choose all possible j ∈ [i, j].

• For each j , choose all possible helices H 1 (i, j , k, l).

• For each H 1 (i, j , k, l), choose all possible H 2 (i , j, k , l ) such that k < i < k < l and l > j .

The compatibility of H 1 (i, j , k, l) and H 2 (i , j, k , l ), ie. whether k < i < k < l and l > j , can be checked directly in the table n×(2∆+1). The partition function

Z K of a single K H ( H 1 , H 2 ) is Z K = e -E K H k B T × Z H 1 × Z H 2
where E K is the energy penalty for constructing a pseudoknot K H . Notice that this equation does not count in the contributions inside the pseudoknot. The partition function Z K of K (i, j) is then

Z K (i, j) = K( H 1 , H 2 )∈K (i,j) Z K .
Again, this partition function counts the pseudoknots and pseudoknots only.

The time complexity of this algorithm is O(n 6 × ∆ 2 ), given by the selection of two helices. The memory complexity is O(n 3 × ∆) since for each K (i, j) we need only to remember its partition function which is then traced back.

Pseudoknots and partially ordered sets

Here we shortly discuss the usability of the partially ordered sets, shortly posets, with respect to the pseudoknots. The idea is to organize the imperfect branches from H ∆ i,j in posets.

Definition 5.2.5 (Poset of helices):

The poset of helices is a graph G(V, E) where:

• V are the imperfect branches H(i, j, k, l);

• E are directed edges from H P (i, j, k, l) to H C (i, j, k , l ) where k ≤ k < l ≤ l
and the k , resp. l is minimal, resp. maximal.

The minimality, resp. maximality of k , resp l implies that H C (i, j, k , l ) should be as similar to H P (i, j, k, l) as possible. Other helices that include H P (i, j, k, l) can be accessed by transitivity.

Example 5.2.2. The example is provided on Figure 5.3. As shown, it does not have to be a tree, as a single child branch H C (i, j, k , l ) may have multiple parents H P (i, j, k, l).

i j i i + 1 j -1 j i i + 2 j -1 j i i + 1 j -2 j H a (i, j, i, j) H b (i, j, i + 1, j -1)
H c (i, j, i + 1, j -1)

H d (i, j, i + 1, j -2)
H e (i, j, i + 2, j -2) The reasoning behind organizing imperfect helices into posets is the automatic incompatibility of H C (i, j, k , l ) with some H 1 (i 1 , j 1 , k 1 , l 1 ) if its parent H P (i, j, k, l)

i i + 2 j -2 j
is incompatible when constructing a pseudoknot. In other words, any of children of H P (i, j, k, l) are not compatible with H 1 (i 1 , j 1 , k 1 , l 1 ) if H P (i, j, k, l) itself is incompatible with it. Since H C (i, j, k , l ) consumes more nucleotides than H P (i, j, k, l), it means that it would produce bigger overlap with H 1 (i 1 , j 1 , k 1 , l 1 )

when constructing a pseudoknot. The opposite is not true, it is still necessary to check the specific branch H C (i, j, k , l ) even if all its parents are compatible with H 1 (i 1 , j 1 , k 1 , l 1 ). However, an incompatible parent H P (i, j, k, l) allows to automatically exclude all H C (i, j, k , l ) from constructing a pseudoknot K H ( H 1 , H C ), potentially reducing a time complexity.

Unfortunately, in practice it is simpler to just allocate a n × (∆ + 1) matrix, since this allows us to directly search only for compatible helices once H 1 (i 1 , j 1 , k 1 , l 1 ) is chosen. The preliminary performed on short sequences (order of 30 nt) with ∆ = 4 have shown three times longer execution times for poset implementation than for the simple matrix. This is mainly due to memory allocation overhead related to the posets aggravated by the fact that a separate poset is needed for each (i, j). Therefore, we decided to use the simple table implementation, which is also discussed in the rest of this Section.

DP scheme for secondary structures with pseudoknots

Once we designed and algorithm to enumerate all pseudoknots K H ( H 1 , H 2 ) with restriction ∆ and to compute the associated partition function, we must establish a DP scheme that allows for sampling of all secondary structures having such pseudoknots. The easiest way to do this is to extend some existing DP scheme (Q, q root , ρ) so it includes the specified pseudoknots. This can be achieved by adding a dedicated table that stores the partition function, and whose values are precomputed by the algorithm that enumerates all pseudoknots.

Since we need to sample the secondary structures in a Boltzmann distribution, we choose the DP scheme employed in McCaskill algorithm as our starting point. We consider a new matrix K, that stores the partition functions Z K i,j of the all pseudoknots spanning over w[i, j] and includes the internal sequences it delimits. That means that for K H ( H 1 , H 2 ), with H 1 (i, j , k, l) and H 1 (i , j, k , l ), we we need to take into an account what happens for the substrings w[k, i ], w[k , l] and w[j , l ]. The values of these sequences internal to pseudoknot depend on how we want the pseudoknot to be scored.

Since each pseudoknot K H (h 1 , h 2 ) can be seen as a special case of a multibranch loop, we might want to evaluate the internal segments as a special case of an internal loop, evaluating for example the structures of a segment w[k, i ] by Z M k,i . However, in pseudoknot the entire w[k, i ] can be unpaired, while Z M k,i does not account for such a case.

To solve this problem, we introduce the new matrix I such that

Z I i,j = Z M i,j + c × (j -i + 1)
where c × (j -i + 1) is the penalty for leaving w[i, j] unpaired in multibranch loop and Z M i,j accounts for all other possibilities. The pknotsRG and pKiss use the default value of 9kcal.mol -1 as a penalty for creating a H-type pseudoknots; we used the same value.

Since K H (h 1 , h 2 ) is composed from h 1 (i, j , k, l) and h 1 (i , j, k , l ), we also need to include the matrix that allows us to determine this helices from H 1 (i, j , k, l) and H 1 (i , j, k , l ) respectively. Since from K i,j we know we target H(i, j, k, l) we have

Z H i,j,k,l such that (k, l, Z H i,j,k,l ) ∈ H δ i,j
and the recursive formula to construct all helices H(i, j, k, l) from H(i , j , k, l) is directly given by Equation 5.1:

Z H i,j,k,l = (k,l,Z H i ,j ,k,l )∈H δ i ,j e E ILG (i,j,i ,j ) k B T × Z H i ,j ,k,l .
We can now define Z K i,j . Since K H ( H 1 , H 2 ), with H 1 (i, j , k, l) = H 1 and H 2 (i, j , k, l) = H 2 , must belong to K (i, j), using Definition 5.2.2 we may specify Z K i,j as:

Z K i,j = k<i <k <l<j <l i<k l <j e -E K k B T × Z H i ,j,k ,l × Z H i ,j,k ,l × Z I k,i × Z I k ,l × Z I j ,l
where E K is the energy penalty for creating K H ( H 1 , H 2 ) and Z H is the partition function of H(i, j, k, l).

The last step is to establish a cases where a pseudoknot may appear. In the case of McCaskill algorithm this is rather simple, as we consider that pseudoknot appears in all cases where a regular base pair does. This means that we can simply duplicate all the rules where Z C i,j appears, then replace it with Z K i,j .

However, it is necessary to proceed with caution in the case of the pseudoknot appearing inside a base pair alone. Since we decided to count its energy as for multibranch loop, we must have to count in the contributions for unpaired stretches w[i + 1, k -1] and w[l + 1, j -1] as in the case of multibranch loop, as well as the penalties a and 2b for opening the multibranch loop itself and for creating two branches within it respectively. If the pseudoknot appears in a multibranch loop with other loops, only the contribution 2b for creating two branches is added.

The complete DP scheme, the modified McCaskill algorithm that includes the pseudoknots, is given below:

Z F i,j = Z F i,j-1 + i≤k≤j-θ-1 Z F i,k-1 × Z C k,j + i≤k≤j-θ-1 Z F i,k-1 × Z K k,j Z C i,j = e -E H (i,j) k B T + i<k<l<j k≤l-θ-1 e -E ILG (i,j,k,l) k B T × Z C k,l + + i<k<j e -a-b k B T Z M i+1,k-1 × Z M1 k,j-1 + i<k<l<j k≤l-θ-1 e -a-2b-c×(j-l+k-i-2) k B T Z K k,l Z M i,j = i<k<j e -c×(k-1) k B T Z M1 k,j + i<k<j Z M i,k-1 × Z M1 k,j Z M1 i,j = e -c k B T × Z M1 i,j-1 + e -b k B T × Z C i,j + e -2b k B T × Z K i,j Z K i,j = k<i <k <l<j <l i<k l <j e -E K k B T × Z H i ,j,k ,l × Z H i ,j,k ,l × Z I k,i × Z I k ,l × Z I j ,l Z I i,j = c × (j -i + 1) + Z M i,j Z H i,j,k,l = (k,l,Z H i ,j ,k,l )∈H δ i ,j e E ILG (i,j,i ,j ) k B T × Z H i ,j ,k,l (5.2)
This DP scheme is depicted on Figure 5.4. As seen, the appearance of the pseudoknot K i,j copies that of a base pair C i,j . To treat pseudoknots in an unambiguous manner we need three supplementary matrices, bringing their total number to 7.

The complexity of this algorithm is largely determined by the complexity of of pseudoknot search, since its complexity is O(n 6 × ∆ 2 ), when compared to the complexity O(n 3 ) of the rest of the algorithm if the length of internal loop unpaired stretches is limited. Likewise, the memory complexity is O(n 3 × ∆) due to storage of partition functions for helices.

F i,j

i j i j-1 j i k-1 k j i k-1 k j C i,j i j i j i k l j i i+1 k-1k j-1 j i k l j M i,j i j i k-1 k j i k-1 k j M 1 i,j i j i j-1 j i j i j K i,j i j i k i' k' l j' l' j I i,j i j i j i j H i,j,k,l i k l j i i' k l j' j
Figure 5.4: The DP scheme for the sampling of the secondary structures including pseudoknots. K i,j and I i,j indicate a substring w[i, j] holding a pseudoknot and a substring internal to a pseudoknot respectively. H i,j,k,l stores imperfect helices h(i, j, k, l). See Section 5.2.4 for further details.

Formalizing the extension

We can now formalize the DP scheme for pseudoknotted structures described in this Section. We denote it by (Q MCv , q C 1,n , ρ MCv ), with MCv standing for Mc-Caskill Variant. The parting point for establishing ρ MCv is ρ MCv is described by Equation 2.20 describing VZualgorithm, which has the same formalization as Mc-Caskill algorithm.

Taking into account how the energy of the proposed pseudoknots is computed, the DP scheme (Q MCb , q C 1,n , ρ MCv ) introduces six new elements:

• A type-H pseudoknot appearing at the level of external loop;

• A type-H pseudoknot appearing inside base pair alone;

• A type-H pseudoknot appearing within a multibranch loop;

• A decomposition of a pseudoknot into its elements, namely imperfect branches H 1 and H 2 and the intermediate segments w[k, i ], w[k , l] and w[j , l ].

• An unpaired segment appearing inside a pseudoknot, not counting imperfect helices;

• A loops appearing inside a pseudoknot, not counting imperfect helices;

This translates into establishing a new constructors for each element:

• λ pk F i,j,k a pseudoknot in an external loop; • λ pk C i,j,k,l a pseudoknot in a base pair; • λ pk M i,j a pseudoknot as a branch in multibranch loop;

• λ pk H i,j,k,l,i ,j ,k ,l a decomposition of a pseudoknot to imperfect helices h 1 (i, j , k, l) ∈ h 1 (i, j , k, l), h 2 (i , j, k , l ) ∈ h 2 (i , j, k , l ). Creating a helix is equivalent to creating a stack/bulge/internal loop -in that case we can reuse the constructor λ ILG i,j,k,l .

Therefore, ρ MCv may be given as

ρ MCv (q F i,j ) =                      ({q F i,j-1 }, λ unp i,j , 1) i≤k<j (k,j)∈S(i,j) ({q F i,k-1 , q C k+1,j-1 }, λ pair i,j,k , 1) i≤k<j (k,j)∈S(i,j) ({q F i,k-1 , q K k+1,j-1 }, λ pk F i,j,k , 1) ρ MCv (q C i,j ) =                                (∅, λ H i,j , E H (i, j)) i<k<l<j (k,l)∈S(i,j) ({q C k,l }, λ ILG i,j,k,l , E ILG (i, j, k, l)) i<k<l<j (k,l)∈S(i,j) ({q K k,l , λ pk C i,j,k,l }, a + 2b + c × (j -l + k -i -2)) i<k<j ({q M i+1,k-1 , q M1 k,j-1 }, λ ML i,k,j , a + b) ρ MCv (q M i,j ) =          i<k<j ({q M1 k,j }, λ ML1hel i,j,k , c × (k -i)) i<k<j ({q M i,k-1 , q M1 k,j }, λ ML2hel i,j,k , 1) ρ
MCv (q M1 i,j ) =            ({q M1 i,j-1 }, λ unp i,j , c) ({q C i,j }, λ p i,j , b) ({q K i,j }, λ pk M i,j , -2b) ρ MCv (q K i,j ) = i<k l <j k<i <k <l<j <l
({q H i,j ,k,l , q H i ,j,k ,l , q I k,i , q I k ,l , q I j ,l }, λ pkH i,j,k,l,i ,j ,k ,l , E K )

ρ MCv (q I i,j ) =    (∅, λ Ipk i,j , c × (j -i + 1)) ({q M i,j }, λ Mpk i,j , 1) ρ MCv (q H i,j,k,l ) = i<i <k l<j <j ({q H i ,j ,k,l }, E ILG (i, j, i , j ))
The search space and the partition function can be formalized for each state q by exploiting the Definition 2.6.1 and Theorem 2.6.2 respectively. The formalism also validates the Equation 5.2 since these relations stem from the Equation 3.11 that was already proven.

The example of formalizing a DP scheme (Q MCv , q C 1,n , ρ MCv ) shows a real strength of the introduced formalism, since it allows to formalize even a states represented by a disconnected substrings, such as represented by the state q H i,j,k,l . This is possible mainly due to the general definitions of q ∈ Q and λ ∈ D that allows them to take as many parameters as necessary, meaning that it is possible to define internal limits of a helix. Consequently, the introduced formalism can be used to express the DP schemes used in algorithms that sample pseudoknotted structures.

Implementation

The above DP scheme was adapted to an algorithm for searching a suboptimal pseudoknots in a program PKnotFind. The program was written in C. For the base of the McCaskill algorithm, we used the implementation employed in VIEN-NARNA library [START_REF] Lorenz | Vi-ennaRNA Package 2.0[END_REF]. The newest iteration of this library allows for a very handy inclusion of an auxiliary grammar, which is the function we used to implement the rules handling the H-type pseudoknots.

While the current iteration of the library allows to add only one additional rule, this is not a problem as both computation of values I i,j and H i,j can be done before computing K i,j . Likewise for traceback of the sampled solution, these values are accessed only from pseudoknots, meaning their values need to be accessed only its traceback. However, since the auxiliary grammar was not included for the stochastic backtrack at the time, we had to adapt its code for the purpose of pseudoknot sampling.

The preliminary version of the program PKnotFind is available at: https://github.com/JurMich/PknotFind

Results

This section provide a few preliminary results testing notably validating the concept of the algorithm. Unfortunately, due to a lack of time we were not able to provide more substantial results such as benchmarks, which is the task that is left for the future.

Comparison with pKiss

We performed the comparison of the test on random sequence w = GCUACGGUCAUCAUCGUAUAGC of 22 nt with pKiss. For pKiss we used the suboptimal mode with energy band (marked as an absolute deviation) of 140kcal.mol -1 . The strategy employed was the strategy P, not returning kissing hairpins. The minimum hairpin length was set to θ = 3 and we also allowed lone base pairs. For PKnotFind, we defined θ = 3 and ∆ = 2. We used uniform energy model, since our main objective is to compare the search space of both software. We then generated 10 7 samples. The PKnotFind predicts 233032 structures. More importantly, we found all 39 pseudoknotted structures that were predicted by pKiss. This points to a complete coverage of the algorithm and DP scheme it is based on, as we do not omit any important pseudoknotted structures. It also predicts a considerable number of pseudoknotted structures, including those with helix length of 1. It also found the structures that pKiss was not able to find due to being restricted only to perfect helices, such as [{{{{{{{]...}}}}}}..}. . We notice that the number of secondary structures, notably of those with pseudoknots, raises exponentially if i) the lone base pairs are allowed and ii) if conditions on helices are relaxed. This is the case even for low ∆. While PKnotFind is able to cover bigger number of secondary structures than pKiss when it comes to the H-type pseudoknots, this does not always have a positive impact. Notably, the lone base pairs have no stabilizing effect on the secondary structure and the construction of the pseudoknot is accompanied by a heavy free energy penalty, meaning that the construction of the pseudoknot of the helices consisting of lone base pairs is inefficient. For this reason, it would be a good idea to always consider the imperfect helices consisting of stacks of length at least 2 base pairs in the case of PKnotFind.

However, the main advantage of our algorithm is the sampling mode, which allows to sample secondary structures with H-type pseudoknots, since the partition function is computed for all possible pseudoknots spanning over w[i, j]. The energy model of the pseudoknot is based on multibranch loop, but the penalty for creating it still remains to be determined.

Complexity Reduction of the PKnotFind algorithm

The algorithm of the PKnotFind is of complexity O(n 6 × ∆ 2 ), which is too high for longer sequences to be analyzed. For example the sequence AY302558.1 /7298-7418 from RFAM family RF00041, a sequence merely 221 nt long, took 10 min to analyze even with the low value of ∆ = 2. The problem partially results from the presence of the pseudoknots with the lone base pairs, which have low overall contribution to the final partition function but due to their important numbers they significantly contribute to the elevation of the total complexity. This can be however extended to the all pseudoknots with high energy in general. The majority of such pseudoknots is composed of helices that present high energies. Consequently, the solution is to exclude such helices before proceeding to the step of building pseudoknots from them is the most straightforward solution.

To filter the helices, we introduce a parameter we call goodness. Definition 5.4.1 (Goodness): Consider an imperfect helix h(i, j, k, l). The goodness g g g of h(i, j, k, l) is the value:

g = log(Z h ) (j -l + k -i + 2)
where Z h is the Boltzmann factor of h(i, j, k, l).

The goodness essentially evaluates the ratio of value brought by h(i, j, k, l) in the form of log(Z h ) when compared to its cost or number of nucleotides it occupies.

Instead of using a table system described in Section 5.2.1, the helices in H ∆ i,j are stored as lists that are ordered according to their goodness, from highest to lowest. For each set H ∆ i,j , we then keep only B helices with best g. Since the value of log(Z h ) is the highest for the perfect helices, the space of the pseudoknots returned by pKiss is included in that of structures returned by PKnotFind as long as B ≥ 1.

The complexity of an algorithm modified in such way is O(n 4 B 2 ). The search for helices keeps its complexity of O(n 3 ∆ 2 ), but since each H ∆ i,j is now restricted to at most B helices, for each combination of i, j, i , j there is possible to have at most B 2 combinations of helices to create the pseudoknots. Consequently, the memory complexity is also reduced to O(n 2 B) due to need to remember only B helices for each couple (i, j) at most. This modification also removes the problem with the lone base pairs, discussed above, because helices with small Boltzmann factors are less likely to be selected during the filtering process. Indeed when we retest the sequence w = GCUACGGUCAUCAUCGUAUAGC we find only |S n | = 132793 structures for the values of ∆ = 2 and B = 11, demonstrating that even for relatively short sequences, there are many helices that are generated even for small values of ∆ and therefore their filtering can greatly improve the speed of the algorithm especially when applied to longer sequences.

The question this solution introduces is how to specify value of B. This question is even more problematic since the number of helices is higher for longer intervals. The low values of B might reduce the search space too drastically (though the most optimal helices are still kept), while large ones will not improve the execution speed too much and will not prevent the generation of less stable pseudoknots. The ideal solution would be an adaptive parameter that would take the length of the interval into an account. However, further analysis must be performed in order to identify the best method.

Pseudoknots and non-redundant sampling

The implementation of non-redundant sampling principle could not be done due to time constraints. In this section we will briefly discuss on how to implement it into the algorithm using the DP scheme (Q MCv , q C 1,n , ρ MCv ).

The implementation can be done in an equivalen way to the one described in Sections 3.3 and 3.4.1. In fact, the latter describes the implementation in the algorithm used by VIENNARNA library, the same algorithm that was used as a starting point for the algorithm described in this Section. The only difference is the complexity of the implementation, since we need to ensure for pseudoknots that they are unique as well.

To ensure this, we proceed in the same manner as for other cases. The nodes of the data-structure d(w, Q MCv , q C 1,n , ρ MCv ) store the partition function of all structures accessible from incomplete structure Z F S I represented by a parse tree T I (w, Q MCv , Q NMCv , q C 1,n , ρ MCv ).

The main difference is the list of new constructors λ that was added along with the associated derivations, which means the extra cases to be handled and consequently increased complexity.

The formalism is useful here since it shows us exactly what new cases were added and must be treated. Likewise any scheme that can be formalized that way is also susceptible to non-redundant sampling, which was globally proven in Section 3.2.1. Therefore the implementation of non-redundant sampling to an algorithm based on DP scheme (Q MCb , q C 1,n , ρ MCv ) is definitely possible and the only constraint is to handle all necessary derivations, which are more numerous due to a presence of extra elements.

Conclusion -Extension on Pseudoknots

In this section we presented the algorithm of the PKnotFind, an extension of an algorithm used in pknotsRG, resp. pKiss. This algorithm enumerates all imperfect helices that do not present longer unpaired stretches and bigger differences between the length of both sides of a helix than ∆ and no bigger deviation between consecutive base pairs is present. These helices are then combined into H-type pseudoknots if they do not overlap.

The complete algorithm to sample stochastically the pseudoknotted structures is generated by incorporating the said algorithm into McCaskill algorithm based on DP scheme (Q MC , q C 1,n , ρ MC ). This is done by enumerating all cases where a pseudoknot can appear, which are identical to the cases of appearance of a base pair. However, the energy constants are different, since in this algorithm, we treat a pseudoknot as a multibranch loop. We therefore needed to duplicate the rules for base pair construction while adjusting the mentioned values.

The practical tests are only preliminary, but they confirmed the expected behavior of the algorithm for T → +∞ with free energy values unaffected, namely the fact Z +∞ = |S n | and the samples being uniformly distributed. The former could be verified for few sequences specifically designed for this purpose, while the latter was also tested for sequences with |S n | < 2.5 × 10 4 by extensive sampling. We also confirm that PKnotFind is able to find the pseudoknotted secondary structures predicted by pKiss, pointing to the correct conceptual scheme behind the software. In addition it is able more pseudoknots than PKnotFind due to the relaxed conditions on the helix. However, not all of these structures seem to be useful, particularly the imperfect helices should contain only the stacks at least 2 base pairs long. This is addressed by a filtering method that removes the helices with higher energies before the step of assembling the pseudoknots. This modification improves time and space complexity of the algorithm and removes less interesting pseudoknots from the search space. However, further analysis on the filtering parameter must be performed, and more tests are necessary to completely validate the PKnotFind software along with the eventual benchmark tests.

The following task that must be completed is the implementation of non-redundant principle into PKnotFind. This is simplified by the formalization of DP-scheme (Q MCv , q C 1,n , ρ MCv ), since it clearly shows newly added derivations and corresponding constructors. This allows to adapt simply the non-redundant version of McCaskill algorithm in VIENNARNA to pseudoknots and of course to adapt the pseudoknot constructing algorithm itself by expanding the non-redundant layer to data-structure d(w, Q MCv , q C 1,n , ρ MCv ). This can be done, seeing that the given DP scheme can be formalized, and therefore it constitutes a perspective into the future.

163

Chapter 6

Conclusion and Perspectives

Our main objective concentrated on the analysis and application of the non-redundant sampling in the domain of the structural RNA bioinformatics. Since in the most of the cases, the redundancy is non-informative, it only hinders the sampling process by slowing it down. We studied the methods of the removal of such redundancy without introducing bias and its application in the domain of structural bioinformatics. We notably studied the efficiency of such approach on the folding landscape modeling and by the comparison with other, classic methods. We also studied the estimation of statistical properties from non-redundant sampling sets and the possible extension of non-redundant sampling on pseudoknotted structures. This study opens the door towards wider application of the non-redundant sampling principle on DP-based problems, both within and outside the field of RNA structural bioinformatics.

We began by establishing the concept of non-redundant sampling and the elements that are related to it. It was demonstrated that for the non-redundant sampling without bias to be possible, it is necessary to remember not only the secondary structures that were sampled, but also the process of their generation. For this reason, we proposed a specific data-structure that stores and handles these informations. This data-structure must be as efficient as possible, since its efficiency greatly determines the overall performance of the non-redundant sampling. We also designed this structure in a manner that it constitutes a separate layer with a limited number of interactions with the algorithm that is extended by a non-redundant sampling. This allows an easy implementation of the principle to any existing algorithms, and when necessary its variables can be adapted for arbitrary precision arithmetics via MPFR library [START_REF] Fousse | MPFR: a multiple-precision binary floating-point library with correct rounding[END_REF].

As the first test, we implemented the non-redundant principle in the algorithm of Saffarian [START_REF] Saffarian | RNA locally optimal secondary structures[END_REF], a combinatorial algorithm that samples saturated secondary structures that are assumed to be locally optimal. Since such structures have a considerable impact on the RNA kinetics due to acting as the kinetic traps, such algorithm is a prime candidate for RNA folding landscape analysis. The generated results demonstrate that the Saffarian algorithm extended by non-redundant sampling samples unique secondary structure more efficiently than the competitors. In addition, the returned models of RNA landscapes were, according to our criteria, closer to the reality than the results obtained by classic sampling methods.

To give users an easy access to the non-redundant sampling principle, we implemented it into the popular ViennaRNA library [START_REF] Lorenz | Vi-ennaRNA Package 2.0[END_REF]. This library employs the DP scheme of McCaskill algorithm. We compared the efficiency of both nonredundant and classical versions of the sampling implemented in the library. With the efficient implementation the non-redundant sampling takes only 10-20% more time than the classic version to generate the same amount of samples. When it comes to unique samples, our algorithm clearly performs better for shorter sequences. For longer sequences, it performs worse at first due to small probability of generating a duplicates and the proportional overhead. However, since the redundant and non-redundant samplings are of the same complexity, it is expected to become more efficient after a sufficient number of the samples is generated.

One may ask whether the redundancy is really devoid of any information, and may point to the problem of the estimation of specific quantities of secondary structures. It is impossible to apply the naive estimator -computing wighted average of all samples -on non-redundant sampling sets. However, we researched an unbiased non-trivial estimator that allows perform such task. We tested the performance of such estimator on some model cases, with demonstrating that our estimator performs better on average that the simple case. However, the problem is the non-trivial estimator does not converge for the entire sampling run. This results from the limited number of unique samples and can be remedied by supposing that the values of estimated statistical quantities of samples whose order exceeds the maximum number of unique structures is zero, making the estimator convergent.

Finally, we tackled the problem of the sampling of pseudoknotted secondary structures. The pseudoknots are usually omitted despite their biological rele-vance due to the impact on the complexity if one searches for pseudoknots [START_REF] Akutsu | Dynamic programming algorithms for RNA secondary structure prediction with pseudoknots[END_REF]. The algorithm used in pknotsRG [START_REF] Reeder | pknotsRG: RNA pseudoknot folding including near-optimal structures and sliding windows[END_REF] and pKiss [START_REF] Janssen | The RNA shapes studio[END_REF] is an efficient algorithm that allows to find the secondary structures with H-type pseudoknots composed of perfect helices. We extended this principle on imperfect helices presenting unpaired bases. The McCaskill algorithm serves as the starting point for our approach to sample secondary structures presenting such pseudoknots. While the algorithm shows it can sample the same pseudoknotted structures as pknotsRG and many more possessing imperfect helices, their number is too high and some of them may be uninformative. This is the specifically the case of the pseudoknots composed from single base pair helices. For this reason we included the filtering of the helices based on the number of nucleotides they occupy and their energy. Such filtering improves space and time complexity of the algorithm and removes pseudoknots containing higher-energy helices from the search space. However, the filtering parameters require further study to most efficiently filter said pseudoknots. Another step that remains to be done is the implementation of the nonredundant sampling principle to the algorithm described above.

To demonstrate the different properties, we also introduced the general language to formalize the DP schemes employed algorithms used in structural RNA bioinformatics. This language is powerful due to its relative simplicity and flexibility that allows to formalize numerous DP schemes with it and prove their properties all at once. As shown, it can be used even to express even the DP schemes employed in algorithms for sampling pseudoknotted structures. It can also efficiently describe the most frequently employed scoring strategies, such as minimization of the energy and partition function limitation.

The final implementations include the software RNANR, written in C with total of 4000 lines, done in collaboration with H. Touzet, that includes non-redundant sampling as well as exhaustive enumeration of the saturated secondary structures, and the choice between double and arbitrary numerical precision arithmetic. The implementation of the non-redundant sampling principle in VIEN-NARNA package was written in C in 1700 lines along with the version including arbitrary precision arithmetics using MPFR library [START_REF] Fousse | MPFR: a multiple-precision binary floating-point library with correct rounding[END_REF]. The relatively low number of lines of code demonstrates the easiness of its implementation particularly well.

The final program, PKnotFind, that is used for sampling of the secondary structures with H-type pseudoknots, was also written in C, uses the algorithm and part of code of VIENNARNA as its base, and has 2000 lines of codes (more for the version containing posets that was ultimately abandoned). Besides these soft-ware, the scripts of total length around 22000 lines of code were created, mostly in Python an R to perform the tests and visualize the results.

Perspective-wise, the non-redundant sampling and the benefits it brings have many possibilities of application in the field of structural RNA bioinformatics. In general, there are many problems where it is necessary to predict a sequence, structure or a composition. Such prediction problems frequently require an access to the suboptimal structures. As long as such problem can be solved using dynamic programming,these problems might be subjected to the the non-redundant sampling principle. It is also useful in the cases where it is necessary to establish a model of the search space, such as we demonstrated in the case of the folding landscapes. This is not restricted to the domain of the secondary structures; for example, the non-redundant sampling might be useful to sample the space of sequences for the purpose of RNA design.

To study the search spaces and similar problems, another class of approaches called parametric inference models [START_REF] Pachter | Parametric inference for biological sequence analysis[END_REF] exists. It represents each biological problem as a graphical model defined by a certain number of parameters. These graphical models are then solved to access the value of these parameters or corresponding space delimited by the range of their values. At first look it may seem that such methods are in competition with non-redundant sampling by basically specifying it by a set of parameters. However, even parametric inference models can benefit from non-redundant sampling. More often than not the parameters feature a certain degree of liberty. Sometimes the number of free parameters can be high, equating to important size of combinations. In this case the sampling can used to extract interesting combinations of parameters, and non-redundant sampling would allow access faster to higher number of them.

But the principle applied for the non-redundant sampling has much richer application than the non-redundant sampling itself. Since the non-redundant sampling principle simply sees the sampled structures as a forbidden structures, it can be easily adapted to not sample certain structures from the beginning, all without introducing the bias between all other structures. This means we can easily sample only certain secondary structures that do not present a certain local substructure, such as specific forbidden loop, or just a specific structure that was blacklisted beforehand, making it potentially a very powerful tool.

The non-redundant sampling has also the potential application outside the field of RNA bioinformatics due to its general concept [START_REF] Andrew | Non-redundant random generation algorithms for weighted context-free grammars[END_REF]. One possible application would concern the simulations of the networks where each branch has maximum capacity. The non redundant sampling would restrict the number of times that a specific path can be taken depending on these capacities, and would perform the simulation of the behavior of such network. This might potentially have a significant use in the domain of the medicine, where different circulatory systems could be simulated this way. This could be helpful to study the circulatory system diseases via the simulations, reducing or removing the necessity of study in humans. It could potentially have the application in other domains as well, such as the urban engineering and traffic networks.

All in all, the non redundant sampling provides a method for optimization of algorithm performance, within the domain of structural RNA bioinformatics or outside of it, and might be the step towards the analysis of the longer RNA sequences. Abstract : Sampling methods are central to many algorithmic methods in structural RNA bioinformatics, where they are routinely used to identify important structural models, provide summarized pictures of the folding landscapes, or approximate quantities of interest at the thermodynamic equilibrium. In all of these examples, redundancy within sampled sets is uninformative and computationally wasteful, limiting the scope of application of existing methods. In this thesis, we introduce the concept of non-redundant sampling, and explore its applications and consequences in RNA bioinformatics. We begin by formally introducing the concept of non-redundant sampling and demonstrate that any algorithm sampling in Boltzmann distribution can be modified into non-redundant variant. Its implementation relies on a specific data structure and a modification of the stochastic backtrack to return the set of unique structures, with the same complexity. We then show a practical example by implementing the non-redundant principle into a combinatorial algorithm that samples locally optimal structures. We use this tool to study the RNA kinetics by modeling the folding landscapes generated from sets of locally optimal structures. These structures act as kinetic traps, influencing the outcome of the RNA kinetics, thus making their presence crucial. Empirical results show that the landscapes generated from the non-redundant samples are closer to the reality than those obtained by classic approaches. We follow by addressing the problem of the efficient computation of the statistical estimates from non-redundant sampling sets. The absence of redundancy means that the naive estimator, obtained by averaging quantities ob-served in a sample, is erroneous. However we establish a non-trivial unbiased estimator specific to a set of unique Boltzmann distributed secondary structures. We show that the non-redundant sampling estimator performs better than the naive counterpart in most cases, specifically where most of the search space is covered by the sampling. Finally, we introduce a sampling algorithm, along with its non-redundant counterpart, for secondary structures featuring simple-type pseudoknots. Pseudoknots are typically omitted due to complexity reasons, yet many of them have biological relevance. We begin by proposing a dynamic programming scheme that allows to enumerate all recursive pseudoknots consisting of two crossing helices, possibly containing unpaired bases. This scheme generalizes the one proposed by Reeders and Giegerich, chosen for its low time and space complexities. We then explain how to adapt this decomposition into a statistical sampling algorithm for simple pseudoknots. We then present preliminary results, and discuss about extensions of the non-redundant principle in this context. The work presented in this thesis not only opens the door towards kinetics analysis for longer RNA sequences, but also more detailed structural analysis of RNAs in general. Non-redundant sampling can be applied to analyze search spaces for combinatorial problems amenable to statistical sampling, including virtually any problem solved by dynamic programming. Nonredundant sampling principles are robust and typically easy to implement, as demonstrated by the inclusion of non-redundant sampling in recent versions of the popular Vienna package.

Universit é Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

1. 1 A

 1 planar representation of chemical structure of double-stranded DNA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 A planar representation of chemical structure of single-stranded RNA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 The five nucleotides occuring within nucleic acids . . . . . . . . . . 1.4 The three levels of RNA structure . . . . . . . . . . . . . . . . . . . . 2.1 The main base pairs occuring within the RNA. . . . . . . . . . . . . 2.2 An example of RNA secondary structure and its different elements shown using different representations. . . . . . . . . . . . . . . . . . 2.3 The decomposition of secondary structures used in Smith and Waterman algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 A loop decomposition of an RNA secondary structure. . . . . . . . 2.5 The principle of Zuker Algorithm . . . . . . . . . . . . . . . . . . . . 2.6 A variant of Zuker (VZu) algorithm used in VIENNARNA library . 2.7 The examples of a folding landscapes . . . . . . . . . . . . . . . . . 3.1 The decomposition employed by Saffarian algorithm . . . . . . . . 3.2 The parameters defined within the Algorithm of Saffarian . . . . . 3.3 The statistics on the secondary structures from the database RNA STRAND . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Decomposition tree T (w, Q ex , q root , ρ ex ). . . . . . . . . . . . . . . . . 3.5 Illustration of the incomplete structure. . . . . . . . . . . . . . . . . 3.6 Simplification of a decomposition tree T (w, Q, q root , ρ) into a datastructure d(w, Q, q root , ρ). . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Linking the parent and children nodes in the data-structure for non-redundant sampling in tree d(w, Q Zuk , q F 1,9 , ρ Zuk ) . . . . . . . . 3.8 The effect of the sorting of list of available flat structures . . . . . . 3.9 The theoretical speed-up for the sampling of two sequences . . . . 3.10 The comparison of the speed gain of RNANR with classical sampling approaches. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 The results of the quality analysis by the different software. . . . . . XV 3.12 An example of folding landscapes between RNANR and RNASLOpt. 3.13 The performance of the non-redundant sampling implementing the hash table organization of the non-redundant layer compared to the usual sampling approach by RNAsubopt . . . . . . . . . . . . 3.14 The performance of the non-redundant sampling implementing the linked list organization of the non-redundant layer compared to the usual sampling approach by RNAsubopt . . . . . . . . . . . . 3.15 The performance of the non-redundant sampling implementing the linked list organization of the non-redundant layer and storing the values using MPFR compared to the usual sampling approach by RNAsubopt. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.16 The performance of the non-redundant sampling implementing the linked list organization compared to the usual sampling approach by RNAsubopt, counted for the number of unique structures.116 4.1 The comparison of efficiency e R and e NR for the estimation of a dotplot matrix D for 10 3 samples. . . . . . . . . . . . . . . . . . . . . 4.2 The histogram of e R -e NR for the estimations of a dotplot matrix D. 4.3 The performance of the estimators D(X ) and D(U) for the two sequences w from the family RF00001 and RF01071. . . . . . . . . . 4.4 The comparison of distributions of distances D i,j for sampling sets X and U generated by classical and non-redundant sampling respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 The evolution of e R and e NR for the number of the samples of X and U of sequence M30199.1/68-167. . . . . . . . . . . . . . . . . 4.6 The distribution of the values B 1 and B 2 obtained for the sequence M30199.1/68-167. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 The evolution of the estimators F(X ) and F(U) for the number of the samples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 An example of a canonical simple recursive pseudoknots represented radially and linearly. . . . . . . . . . . . . . . . . . . . . . . . 5.2 An example of imperfect helix contained by h(i, j, i + 8, j -8). . . . 5.3 The organization of helices in poset. . . . . . . . . . . . . . . . . . . 5.4 The DP scheme for the sampling of the secondary structures including pseudoknots. . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: A planar representation of chemical structure of double-stranded DNA. Each strand is a chain of nucleotides, connected by phosophodiester bond and composed of Phosphate (purple), 2-Deoxyribose (black) and a nitrogenous base. The base has four main variants -Adenine (A, green), Guanine (G, Yellow), Cytosine (C, Blue) and Thymine (T, brown). The two strands are held together by hydrogen bonds between nitrogenous bases following the base-pairing rules. The 5' end and 3' end are determined by the orientation of the 2-Deoxyribose.

Figure 1 . 2 :

 12 Figure 1.2: A planar representation of chemical structure of single-stranded RNA.In this case, the sugar (black) is ribose and Uracil (U, orange) substitutes Thymine (T, brown).

Figure 1 . 3 :

 13 Figure 1.3: The five nucleotides occuring within nucleic acids. Thymine does not occur within RNA while Uracil is specific to it.

Figure 1 . 4 :

 14 Figure 1.4: The three levels of RNA structure. The primary structure (A) is essentially a chain of nucleotides. The secondary structure (B) also includes the hydrogen bonds between the nucleotides according to valid base pairing rules. The tertiary structure (C) takes into account all intramolecular forces. Molecule reference: PDB ID 1FIR [13, 104, 8].The secondary structure taken from RNA STRAND[START_REF] Andronescu | RNA STRAND: the RNA secondary structure and statistical analysis database[END_REF] and visualized by VARNA[START_REF] Darty | VARNA: Interactive drawing and editing of the RNA secondary structure[END_REF]. The tertiary structure created by PYMOL[START_REF] Schrödinger | The PyMOL Molecular Graphics System[END_REF].

Definition 2 . 1 . 1 (

 211 Sequence): An RNA sequence w w w is a chain of nucleotides of length n over the alphabet Ω = {A, C, G, U}.

Definition 2 . 1 . 5 (

 215 Secondary structure): An RNA secondary structure, or simply secondary structure S(i, j) is a set of base pairs (a, b) (a, b) (a, b) within a substring w[i, j] (∀a, b, i ≤ a < b ≤ j). We have, ∀i, j, S(i, j) ⊆ φ pair (w). The secondary structure of entire sequence w w w is denoted by S(1, n) = S.

Figure 2 . 1 :

 21 Figure 2.1: The main base pairs occuring within the RNA. These base pairs consist of hydrogen bonds (red dashed lines) between the atoms of nucleotides. The Adenosine -Uridine (left) and Guanosine -Cytidine(center) base pairs are also called Watson-Crick.The Guanosine-Uridine (right) base pair is also called wobble and is less frequent due to slightly lower energy.

Figure 2 . 2 :

 22 Figure 2.2: An example of RNA secondary structure and its different elements shown using different representations: radial (A), circular (B), linear (C) and dot-bracket (D).The different elements are an external loop (cyan), hairpin (red), stem/helix (blue), an internal loop (grey), a bulge (white) and a multibranched loop (yellow). Representations A-C were created using VaRNA[START_REF] Darty | VARNA: Interactive drawing and editing of the RNA secondary structure[END_REF].

Figure 2 . 3 :

 23 Figure 2.3: The decomposition of secondary structures used in Smith and Waterman algorithm.The base j can be either unpaired or paired to some k from w[i, j -θ -1]. Both cases produce shorter substrings on which this process can be recursively applied.

Definition 2 . 3 . 2 (

 232 Suboptimal structure): Suppose each secondary structure has a certain score S(S) and the objective is to minimize S. We call a suboptimal structure every S T that S(S) > min S T ∈Sn (S(S))

Figure 2 . 4 :

 24 Figure 2.4:A loop decomposition of an RNA secondary structure. The different loops are an external loop (cyan), hairpin (red), stem/helix (blue), an internal loop (grey), a bulge (white) and a multibranched loop (yellow). The secondary structure was created using VaRNA[START_REF] Darty | VARNA: Interactive drawing and editing of the RNA secondary structure[END_REF].

Figure 2 . 6 :

 26 Figure 2.6: A variant of Zuker (VZu) algorithm used in VIENNARNA library. F i,jand C i,j are functionally identical to W i,j , respectively V i,j from the original version of Turner's Algorithm. M i,j and M 1 i,j score the section of a multibranched loop within w[i, j] containing at least one helix and exactly one helix respectively. M i,j is split to a section that contains exactly one helix possibly preceeded by unpaired bases and another that is either completely unpaired or contains at least one another helix. M 1 i,j contains one helix potentially followed by unpaired bases.

Figure 2 . 7 :

 27 Figure 2.7: The examples of a folding landscapes. The folding landscape for Yeast Phenylalanine tRNA (PDB ID 1TRA, left)[START_REF] Westhof | Restrained refinement of the monoclinic form of yeast phenylalanine transfer RNA. Temperature factors and dynamics, coordinated waters, and base-pair propeller twist angles[END_REF] and Class II Amynoacyl tRNA synthetase (1ASY, right)[START_REF] Ruff | Class II aminoacyl transfer RNA synthetases: crystal structure of yeast aspartyl-tRNA synthetase complexed with tRNA (asp)[END_REF]. The first and second reference are respectively MFE structure and a random high-energy structure distant from MFE structure. The picture was generated by an Rscript script that will be possibly available with the future releases of VIEN-

Figure 3 . 1 :

 31 Figure 3.1: The decomposition employed by Saffarian algorithm. The substring w[i, j]is treated depending on whether (i -1, j + 1) is present or not. If it is absent, all flat structures p(i, j) existing for w[i, j] are considered. If it exists, all p(i, j) except H(i, j, a) are considered along with the possibility of extending the helix H(i -a + 1, j + a -1, a) by a single base pair (i, j). The helices are marked by blue.

2 .Figure 3 . 2 :

 232 Figure 3.2:The parameters defined within the Algorithm of Saffarian. ν denotes the maximum base pair span, which allows to define whether we allow the long-rande interactions or not. θ represents the minimum length of a base pair, while π is the maximum length of a hairpin. γ denotes the maximum number of the branches within the multibranch loop.

Figure 3 . 3 :

 33 Figure 3.3: The statistics on the secondary structures from the database RNA STRAND. The statistics are calculated for the classes of sequences separated by their maximum length. (A) The percentage of the structures from RNA STRAND having at maximum the indicated number of branches within each of their multibranch loop. Here this number also include a closing base pair. (B) The proportion of helices that have at least the indicated length. The statistics show that both multibranch loops with more than 6 branches and the lone base pairs are uncommon.

Figure 3 . 5 . 7 }

 357 The sequence w = ACCCAAAGGAGA contains an incomplete structure constituted of base pairs (2, 11) and (3, 9) and an underived state q C 4,8 . We may choose among two possible derivations: creating a stack of base pairs (3, 9) and (4, 8).

8 ∅. 7 ∅

 87 The former will lead to an incomplete structure S I1 = {(2, 11),[START_REF] Arrhenius | Über die dissociationswärme und den einfluss der temperatur auf den dissociationsgrad der elektrolyte[END_REF][START_REF] Biebricher | Structural analysis of self-replicating RNA synthesized by Qβ replicase[END_REF], (4, 8)} and a state q C 5,7 on which only the derivation by q C is possible, giving S 1 = {(2, 11),[START_REF] Arrhenius | Über die dissociationswärme und den einfluss der temperatur auf den dissociationsgrad der elektrolyte[END_REF][START_REF] Biebricher | Structural analysis of self-replicating RNA synthesized by Qβ replicase[END_REF], (4, 8)}. λ H 4,8 leads to a complete structure and S 2 = {(2, 11), (3, 9)}. Therefore

  number of structures N Result: Sset a set of unique secondary structures M DP ← Precompute_DP() # Pre-computes and fills all DP-matricesSset←∅ W ← ∅ d ← N(q, 0, ∅) # constructs node N(id, 0, Nparent) in T(Q, qroot, ρ) ptr ← qroot Z S I ← Z while len(Sset) < N do S ← unfolded while q do λ list ← list_deriv(ptr) Ztemp = 0 r ← random(0,1) # random float between 0 and 1 Nodenext ←null # stores the child of ptr while λ in λ list and Ztemp < r do Nodenext = get_child(λ, ptr, W) # gets child for case λ of ptr, returns null if node does not exist Ztemp = Ztemp+get_probability(M DP , ptr, λ, Nodenext) Z S I ← Z S I ×get_probability(M DP , ptr, λ, FS I q , Nodenext) S ←add_object(λ)# adds element associated to λ List deriv ← get_all_states(ptr, λ) # lists q 1 , • • • , q u for q λ --------→ e -c λ /k B T {q 1 , . . . , q u } if len(List deriv ) > 1 then Push List deriv [2 :] to W if not Nodenext and List deriv = ∅ then Nodenext ←N(λ, 0, ptr) ptr = Nodenext Push S to Sset count = count + 1 ptr ← traceback(ptr)# traceback to root, update nodes

Figure 3 . 9 :

 39 Figure 3.9: The theoretical speed-up for the sampling of two sequences. The theoretical speed-up computed for the sequence of 5S ribosomal RNA of Thermoplasma acidophilum 123 nt long (A) and for Telomerase of RNA of Tetrahymena silvana (154 nt, B). Both sequences were taken from RNA STRAND database. Note the exponential increase with the size of U

  The results are shown on Figure 3.10.

Figure 3 . 10 :

 310 Figure 3.10: The comparison of the speed gain of RNANR with classical sampling approaches.

Results 3 . 3 . 2 .

 332 The results are shown on Figure 3.11. As expected, RNANR most frequently identifies both structures due to being used in their generation (Figure 3.11A

Figure 3 .

 3 Figure 3.11: The results of the quality analysis by the different software. The successful identification shows RNANR was the most successful with retrieval of two metastable structures (A). This is expected due to RNANR being used to create the bistable sequences. However the lowest switching time is returned with better proportion than the retrieval of the structures (B). This hints at statistically better performance of RNANR when it comes to the analysis of RNA landscape. For (B), the 'None' category indicates the case where no software was able to retrieve both structures.

Figure 3 . 12 :

 312 Figure 3.12: An example of folding landscapes between RNANR and RNASLOpt. For the RNA sequence (A) the metastable (B) and MFE (D) locally optimal secondary structures are generated, with helix in red being presend only for SRNANR. The landscapes generated by RNANR and RNASLOpt (C) are similar and consist of two parts, but for

Figure 3 . 13 :

 313 Figure 3.13: The performance of the non-redundant sampling implementing the hash table organization of the non-redundant layer compared to the usual sampling approach by RNAsubopt. The total number of the samples is considered. The four sequences from respective families RF00174 (A), RF00061 (B), RF00001 (C) and RF01071 (D) were sampled for secondary structures by a redundant (blue) and non-redundant (light purple) method. The lines represent the linear approximation of the graphs. The dark purple number indicates the speed ratio between redundant and non-redundant sampling, here indicating the overhead of 50%.

Figure 3 . 14 :

 314 Figure 3.14: The performance of the non-redundant sampling implementing the linked list organization of the non-redundant layer compared to the usual sampling approach by RNAsubopt. The total number of the samples is considered. The four sequences from respective families RF00174 (A), RF00061 (B), RF00001 (C) and RF01071 (D) were sampled for secondary structures by a redundant (blue) and non-redundant (light purple) method. The lines represent the linear approximation of the graphs. The time complexity overhead is in this case reduced to 10%-20%.

Figure 3 . 15 :

 315 Figure 3.15: The performance of the non-redundant sampling implementing the linked list organization of the non-redundant layer and storing the values using MPFR compared to the usual sampling approach by RNAsubopt. The total number of the samples is considered. The two sequences from respective families RF00174 (left) and RF00061 (right) were sampled for secondary structures by a redundant (blue) and non-redundant (light purple) method. Using the arbitrary precision arithmetic considerably raises the time complexity overhead more than doubling the total time complexity.

Figure 3 . 16 :

 316 Figure 3.16: The performance of the non-redundant sampling implementing the linked list organization compared to the usual sampling approach by RNAsubopt, counted for the number of unique structures. Only the unique samples are considered. The four sequences from respective families RF00174 (A), RF00061 (B), RF00001 (C) and RF01071 (D) were sampled for secondary structures by a redundant (yellow) and nonredundant (green) method. The lines represent the linear approximation of the graphs. For the shorter sequences A and C, the non-redundant sampling performs better due to high number of the duplicates generated by the usual sampling approach. For the sequences B and D, it performs worse due to the complexity overhead and little probability to generate duplicates for such long sequences at the beginning.

Definition 4 . 1 . 1 (Example 4 . 1 . 1 .

 411411 Feature function):A feature function F : S n → R that associates to each secondary structure S a certain numerical value. F returning whether a base pair (i, j) is present in the secondary structure S:

Definition 4 . 1 . 3 (

 413 Non-redundant estimator): We define by non-redundant sampling set estimator F(U) the estimator

Example 4 . 1 . 3 .

 413 We return to the example of 100 identical RNA secondary structures with a single base pair and 50 identical molecules with two base pairs, where |U| = 2. The Equation 4.1.2 becomes N BP (U) = 1 2 (N BP ( S 1 ) + N BP ( S 2 ) -N BP ( S 1 ) × p( S 1 | U)) + p( S 1 | U) × N BP ( S 2 ))

Theorem 4 . 1 . 1 (Proof 4 . 1 . 1 .

 411411 Unbiasedness): The non-redundant estimator is unbiased. Consider the conditional expectation E U (F( S)) that E U ( F( S)) = S∈Sn p( S | U) × F( S).

  e R = ||D -D(X )|| n(n -1) and e NR = ||D -D(U)|| n(n -1) where e R and e NR are respectively the redundant and non-redundant error and ||a|| = √ a 2 . The results were classed according to coverage (reminder: Cov(U) = S∈U Z S Z

Results 4 . 2 . 1 .

 421 The results are depicted on Figure4.1. For all Figures, the line

  marks

Results 4 . 2 . 2 .

 422 The results for both results are shown on Figure4.3. For both sequences w we notice the better overall performance of D(U) than that of D(X ) for comparable number of samples and execution time, since D(U) show lower error for the same execution time/number of samples. The results for the number of samples is expected, since the non-redundant sample U holds more information and the associated estimator D(U) shows smaller variance. The better performance of D(U) demonstrates that the non-redundant sampling can be implemented efficiently into the associated DP schemes with low complexity overhead.

Figure 4 . 2 :

 42 Figure 4.2: The histogram of e R -e NR for the estimations of a dotplot matrix D. The red dotted line denotes e R -e NR = 0. The positive values indicate the cases for which D(U) performs better than D(X ). The coloring indicates the coverage intervals Cov of the sampling set of w for which the estimators were calculated.

Figure 4 . 3 :

 43 Figure 4.3: The performance of the estimators D(X ) and D(U) for the two sequences w from the family RF00001 (left) and RF01071 (right). The graphs show the evolution the errors e R (blue) and e NR (orange) in function of the number of samples (top) and execution time (bottom).

Definition 4 . 2 . 2 (

 422 Graph distance):A minimum graph distance d i,j (S) is the minimal distances between two nucleotides w[i] and w[j] observed for some S ∈ S n , with S represented by G S (V, E).

Example 4 . 2 . 1 .

 421 For a secondary structure S = {(1, 10)} of a sequence w with length n = 10, the minimum graph distance is d(2, 9)(S) = 3, following the path (2)-(1)-(11)-[START_REF] Biebricher | Structural analysis of self-replicating RNA synthesized by Qβ replicase[END_REF].By definition, all nucleotides w[i], w[j] neighboring in backbone or in a base pair have a d i,j (S).

Results 4 . 2 . 3 .

 423 We plot the histograms of S n , X and U in function of distance classes, one class constituting a specific distance returned for given sampling set. These plots are shown on Figure 4.4. The plot shows the superposition of the distribution of the distance classes for X and S n or reference (Figure 4.4 left) and U and S n (Figure 4.4 right). Both distributions are very close to the reference. The precise results are shown in Table 4.1. The errors (3 rd column) were computed by: e R = ||D 16,82 (S n ) -D 16,82 (X )|| and e NR = ||D 16,82 (S n ) -D 16,82 (U)||.

Definition 4 . 2 . 3 (

 423 Weighted graph distance matrix):A weighted graph distance matrix M for a sequence w of length n is a matrix n × n containing the values:M i,j = S∈Sn d i,j (S) × Z S Z .

Figure 4 . 4 :

 44 Figure 4.4: The comparison of distributions of distances D i,j for sampling sets X and U generated by classical and non-redundant sampling respectively. The expected distribution is superposed to that for sampling sets X (left) and U (right). Both distributions are close to the expected one.

Results 4 . 2 . 4 .

 424 The results for e R and e NR , respectively B 1 and B 2 , are shown on Figures 4.5, resp.

Figure 4 . 6 .

 46 The comparison of errors e R and e NR indicates that M(U) approaches the reference, M, more quickly than M(U). More importantly, the variation of M(U) is smaller than for M(U).

Figure 4 . 5 :

 45 Figure 4.5: The evolution of e R and e NR for the number of the samples of X and U of sequence M30199.1/68-167. The non-redundant estimator F(U), marked in orange, shows smaller variations than its classic variant F(X ), blue, and in general seems to show being close to the true value.

Figure 4 . 6 :

 46 Figure 4.6: The distribution of the values B 1 and B 2 obtained for the sequence M30199.1/68-167. Blue histogram shows the distribution of number of samples for B 1 (e N R < e R for the first time for given samples X and U), the green B 2 (e N R < e R from this point on). The blue and green dashed lines how both the respective average for B 1 and B 2 .

Example 4 . 2 . 2 .

 422 The structure '((...(...)))' becomes '[[]]'.

Figure 5 . 1 :

 51 Figure 5.1: An example of a canonical simple recursive pseudoknots represented radially (a) and linearly (b). Two perfect helices H 1 (i, j, lh) and H 2 (i , j , lh ) cross each other, creating three internal segments I 1 , I 2 , I 3 .

Figure 5 . 3 :

 53 Figure 5.3: The organization of helices in poset. The imperfect branches where the helices it contains can be encompassed in another are connected by an edge.

  and internal segments w[k, i ], w[k , l] and w[j , l ]. • λ Ipk i,j an unpaired segment in pseudoknot; • λ Mpk i,j a loop in pseudoknot.

Results 5 . 3 . 3 .

 533 The pKiss has returned 11415 structures, 39 of which present Htype pseudoknots with perfect helices. Even with lone base pairs allowed the minimum length of all helices in returned pseudoknots was 2. The pseudoknotted structure having the lowest energy was [[.{{{]]......}}}....., with -1.5kcal.mol -1 .

Titre:

  Echantillonage sans remise en Bioinformatique des Acides RiboNucl éiques Mots cl és : ARN, programmation dynamique, échantillonage non-redonddante, algorithme combinatoire, mod èle thermodynamique, cin étique R ésum é : Un échantillonnage statistique est central à de nombreuses m éthodes algorithmiques pour la bioinformatique structurale des ARNs, o ù ils sont couramment utilis és pour identifier des mod èles structuraux importants, fournir des r ésum és des espaces de repliement ou approcher des quantit és d'int ér êt dans l' équilibre thermodynamique. Dans tous ces exemples, la redondance dans l'ensemble échantillonn é est non-informative et inefficace, limitant la port ée des applications des m éthodes existantes. Dans cette th èse, nous introduisons le concept de l' échantillonnage non-redondante et nous explorons ses applications et cons équences en bioinformatique des ARN. Nous commenc ¸ons par introduire formellement le concept d' échantillonnage non-redondante et nous d émontrons que tout algorithme échantillonnant dans la distribution de Boltzmann peut être modifi é en une version nonredondante. Son impl émentation repose sur une structure de donn ées sp écifique et la modification d'une remont ée stochastique pour fournir l'ensemble des structures uniques, avec la m ême complexit é. Nous montrons alors une exemple pratique en impl émentant le principe d' échantillonnage non-redondant au sein d'un algorithme combinatoire qui échantillonne des structures localement optimales. Nous exploitons cet outil pour étudier la cin étique des ARN, mod élisant des espaces de repliement g én ér és à partir des structures localement optimales. Ces structures agissent comme des pi èges cin étiques, rendant leur prise en compte essentielle pour analyser la dynamique des ARN. Des r ésultats empirique montrent que des espaces de repliement g én ér és à partir des échantillons non-redondants sont plus proches de la r éalit é que ceux obtenus par un échantillonnage classique. Nous consid érons ensuite le probl ème du calcul efficace d'estimateurs statistiques à partir d' échantillons non redondants. L'absence de la redondance signifie que l'estimateur naïf, obtenu en moyennant des quantit és observ és dans l' échantillon, est eronn é. Par contre, nous établissons un estimateur non-trivial non-biais é sp écifique aux échantillons non-redondants suivant la distribution de Boltzmann. Nous montrons que l'estimateur des échantillons non-redondants est plus efficace que l'estimateur naïf, notamment dans les cas o ù la majorit é des l'espace de recherche est échantillonn é. Finalement, nous introduisons l'algorithme d' échantillonnage, avec sa contre-partie non-redondante, pour des structures secondaires pr ésentant des pseudonoeuds de type simple. Des pseudonoeuds sont typiquement omis pour des raisons d'efficacit é, bien que beaucoup d'entre eux poss èdent une grande importance biologique. Nos commenc ¸ons par proposer une sch èma de programmation dynamique qui permet d' énum érer tous les pseudonoeuds compos és de deux h élices pouvant contenir des bases nonappari és qui s'entrecroisent. Ce sch èma g én éralise la proposition de Reeders et Giegerich, choisi pour sa base complexit é temporelle et spatiale. Par la suite, nous expliquons comment adapter cette d écomposition à un algorithme d' échantillonnage statistique pour des pseudonoeuds simples. Finalement, nous pr ésentons des r ésultats pr éliminaires et nous discutons sur l'extension de principe non-redondant dnas ce contexte. Le travail pr ésent é dans cette th èse ouvre non seulement la porte à l'analyse cin étique des s équences d'ARN plus longues, mais aussi l'analyse structurale plus d étaill ée des s équences d'ARN en g én éral. L' échantillonnage non-redondant peut être employ é pour analyser des espaces de recherche pour des probl èmes combinatoires susceptibles à l' échantillonnage statistique, y inclus virtuellement tous probl èmes solvables par la programmation dynamique. Les principes d' échantillonnage non-redondant sont robustes et typiquement faciles à impl émenter, comme d émontr é par l'inclusion d' échantillonnage non-redondant dans les versions r écentes de Vienna package populaire. Title : Non-redundant sampling in RNA bioiformatics Keywords : RNA, dynamic programming, non-redundant sampling, combinatorial algorithms,thermodynamic model, kinetics

An external loop. It

	external loop per sequence. On radial representation, it is the linear base
	(cyan on Figure 2.4). It also contains one or more opening base pairs (k x , l x )
	for following substructures unless w is unfolded.
	• A hairpin. A stretch of unpaired nucleotides delimited by a single closing base pair (i, j). It does not branch any further and consequently it repre-
	sents a terminal point for given branching (red on Figure 2.4). Its energy is
	denoted by E H (i, j).
	• A stem. A set of consecutive base pairs (i, j) and (i + 1, j -1) = (k, l) -a closing base pair is immediately followed by opening base pair. A multiple
	consecutive stems form a helix (blue on Figure 2.4).
	Definition 2.4.1 (Helix): We call by H(i, j, lh) H(i, j, lh) H(i, j, lh) a helix consisting of the base pairs (i, j), (i + 1, j -1) . . . (i + lh -1, j -lh + 1) where lh is the
	length of H(i, j, lh) and i + 2lh -2 + θ < j. There cannot be any unpaired base within w[i, i + lh -1] and w[j -lh + 1, j].
	defines the branches and unpaired nucleotides that are
	not enclosed by any base pair. Consequently, there is no closing base pair
	in an external loop -they are delimited by w[1] and w[n]. There is only one

  The basin B S L B S L B S L of locally optimal structure S L is the section of RNA folding landscape containing all structures S basin S L such that there exists a pathS basin S L → S 1 → . . . → S u → S

	Definition 2.8.1 (Basin): where the energies follow
	E S basin S L

≥ E S 1 ≥ . . . ≥ E Su ≥ E S L

In other words, it is every point that is, on folding landscape, in the direction of the gradient towards S L . Consequently, any secondary structure within B S L should spontaneously transform into S L . Therefore we might consider an approximation of RNA landscape by considering only S L . The basins are connected by the secondary structures represented by saddle points. Definition 2.8.2 (Saddle point): A saddle point S peak S peak S peak between two neighboring basins B S LA and B S LB is a structure S peak such that S peak ∈ B S LA ∩ B S LB and

  However, the basins B S LA and B S LB are more likely to be neighboring the more S LA and S LB are similar. The best course of action is to suppose that the B S LA and B S LB are neighboring if they are similar enough. The locally optimal structures B S LA and B S LB are said to be connected if d(S LA , S LB ) < c d where d is the distance function evaluating the similarity of S LA and S LB and c d is the user defined threshold. The space of all S Ln with B S L connected to B S LA is denoted

Definition 2.8.3 (Connected locally optimal structures):

Table 3 .

 3

1: The comparison whether the saturated secondary structures are locally optimal. The category denoted 'Within search space' marks the cases where S grad ∈ S SAT .

  Zuk } used in Zuker algorithm, the example of a incomplete structure S I for the sequence w=ACCCAAAGGAGA contains two base pairs (2, 11) and (3, 9), and a state q C 4,8 with no λ. Two constructors, λ ILG 2,9,4,8 and λ H 3,9,4,8 , can be attributed to this state to generate two complete secondary structures. S I is constituted from two complete structures: one where 4 th and 8 th nucleotide is paired and one where not.

			4 4 4	8 8 8
			4 4 4	8 8 8
	1 1	12 12 12
	Figure 3.5:	Illustration of the incomplete structure.	Under DP scheme
	{Q Zuk , q F 1,12 , ρ

  1/68-167 from the RF00001 family. The reference is obtained by sampling |X ref | = 10 6 structures. For testing the non-redundant and classical sampling, we generate |X | = |U| = 10 3 samples. We look specifically for the distribution of D 16,82 of nucleotides w[16] and w[82].

Table 4 . 1 :

 41 10 4 samples. For 132 Experiment type Distance classes Distance from D 16,82 (S n ) The comparison of distributions of minimum graph distance d 16,82 for 10 3 samples. Value classes denote the number of different distances returned for given set of samples.

	D 16,82 (S n )	36	0
	D 16,82 (X )	18	0.01813
	D 16,82 (U)	22	0.01496

The results for the equivalent analysis, but the non-redundant sampling being implemented by using the linked lists instead, is depicted on Figure3.14. In this case, the time complexity overhead is reduced to much more manageable 10%-20%. This is caused by two factors. First, the linked list necessitates to allocate

Acknowledgements

Keep in mind however that the above examples are not generalizable and as Figure 4.1 shows, the estimator D(U) can perform worse than its classical counterpart in some cases.

We then followed the evolution of the values of both estimator with regards to the expected value F SHA (S n ). Results 4.2.5. The results of this analysis are depicted on Figure 4.7. For all four sequences, the both estimators approach the expected value after 10 4 samples. It is interesting to observe however that in case of the sequence from RF00001 (Figure 4.7A), F(X ) approaches the value quicker than the non-redundant estimator F(U). In all other cases our estimator approaches to the expected value F SHA (S n ) quicker than F(X ).

Another interesting case is the sequence from the family RF00061 (Figure 4.7C). In this case the value F(U) starts to diverge from the expected value. This can be explained either by the rounding errors, or by the inconsistency of the proposed non-redundant estimator. Since F(U) is not consistent it may diverge once U starts to get too large. Note that in this case j -l = k -i is not necessarily true. The imperfect helices h(i, j, k, l) having identical (i, j) and (k, l) may be brought together in common sets.

A B C D

Definition 5.2.2 (Imperfect branch): An imperfect branch H(i, j, k, l) with i < k < l < j is a set of all helices h(i, j, k, l).

The partition function of H(i, j, k, l) is the sum of Boltzmann factors of h(i, j, k, l):

From H(i, j, k, l) we can establish the definition of H-type pseudoknots.

Definition 5.2.3 (H-type pseudoknot):

By abuse of notation we can write

This definition is very general and encompasses all possible H-type pseudoknots. The canonical simple recursive pseudoknots are their specific subcategory. However, only searching all imperfect helices h(i, j, k, l) would induce at least O(n 4 ) memory complexity to store them all.

To restrict the number of h, we introduce the parameter ∆ that implies:

• For any two consecutive base pairs (i 1 , j 1 ) and (i 2 , j 2 ) such that

• w[i u , j u ] a substring of w[i, k] or w[l, j] that is completely unpaired may be at most of length ∆.

Only the helices that fulfill the above criteria are considered for pseudoknots

Example 5.2.1. Here we suppose an imperfect helix h(i, j, k, l) ∈ H(i, j, k, l) with k = i + 8 and l = j -8, shown on Figure 5.2. This helix is a subset of six base pairs h(i, j, k, l) = {(i, j, (i + 1, j -1), (i + 3, j -5), (i + 4, j -6), (i + 7, j -7), (k, l)}.

If we consider ∆ = 3, then helix in h(i, j, k, l) is considered since there at most 3 consecutive unpaired nucleotides and the maximum skew is between pairs (i + 1, j -1) and (i + 3, j -5), where

On the other hand, if ∆ = 2 then helix in h(i, j, k, l) will be excluded due to an unpaired stretch of length 3 it contains, and h(i, j, k, l) will not be in H(i, j, k, l). ki + 1 = 9 jl + 1 = 9 i i+1 i+3 i+4 i+7 k=i+8 l=j-8 j-7 j-6 j-5 j-1 j H(i, j, k, l) We can verify that this helix fulfills ∆ = 2, and therefore in that case is considered as an element of

Such a restriction is justified since it means that one base pair of h(i, j, k, l) cannot be too distant from the others and prevents its general asymmetry as well. This is important since we are not looking what happens for the unpaired bases between base pairs, meaning too much unpaired base pairs in h(i, j, k, l) would leave too many of them in the final structure.

In the terms of complexity, such a restriction means that instead of having at maximum n 2 possible imperfect helices for each couple (i, j), we have at maximum

Toy Examples

The first test that we performed were destined to verify whether the algorithm covers the H-type pseudoknots K( H 1 , H 2 ) we are interested in. Since the entirety of sequences w of the search space is hard to enumerate for long w, we performed these tests on artificial sequences designed specifically for them. For this reason we assume ∆ = 2 as well.

The easy way to test the completeness of the DP scheme (Q MCb , q C 1,n , ρ MCv ) is to assume a temperature T → +∞ without it affecting the energy values, in which case we have, ∀S:

T →+∞ = 1 which implies uniform distribution and Z +∞ = |S n |. For simpler sequences, it is simple to count manually the number of possible secondary structures, including those with K( H 1 , H 2 ), then validate the results by computing the partition function when T → +∞.

The first sequence for which the test was performed is w = CACAAGAG. The reason we provide the results for this short sequence is it is explicit to enumerate all possible secondary structures that satisfy the given conditions, notably the minimum distance between base pairing nucleotides θ = 3 and ∆ = 2:

These counts point to an uniform distribution of Z +∞ S as stated. Therefore for the simplest examples, the PKnotFind returns expected solutions.

The second interesting sequence we tested is

While the enumeration of search space structures is not as explicit as in previous example, it is still easy to enumerate all secondary structures satisfying the given conditions. In this case |S n | = 33. In this case we sampled N = 10 6 samples.

Results 5.3.2. We found all 33 structures with PKnotFind including the pseudoknotted structures. Notably, we also found the secondary structure [.{..{....]..}}, which would not be found if the we searched only for the perfect helices. This structure also satisfies ∆, with the maximum number of an unpaired bases in helix being 2 as well as the difference between the length of both sides of helices

The results also show that the number of secondary structures raises very quickly if even the simplest pseudoknots are considered and ∆ is small. For example for the random sequence w = GCUACGGUCAUCAUCGUAUAGC we obtain Z +∞ = |S n | = 233032 structures. This makes validation for longer sequences and especially the cases where pseudoknots are combined with multibranch loops difficult, notably to check for uniform distribution since it takes too many samples to have enough samples for each structure to verify. However, we managed to find a few sequences w with |S n | < 2.5 × 10 4 and for N = 10 7 to verify the uniformity of generated distribution as well as that the number of the structures equals to Z +∞ = |S n |, which supports the validity of PKnotFind.

We also benchmarked the preliminary version of PKnotFind used for tests with one based on posets (see Section 5.2.3). The test for sequences of around 30 nt and ∆ = δ have shown thrice as high execution time for the latter version. For this reason we abandoned the development of the version of PKnotFind with posets and stayed with the current one.