<figure><img src="image_1.png" /></figure>

HAL
open science

# Ensemble Algorithms and Analytic Combinatorics in RNA Bioinformatics and Beyond

Yann Ponty

<figure><img src="image_2.png" /></figure>

## To cite this version:

Yann Ponty. Ensemble Algorithms and Analytic Combinatorics in RNA Bioinformatics and Beyond. Bioinformatics [q-bio.QM]. Université Paris-Saclay, 2020. tel-03219977v2

HAL Id: tel-03219977

https://theses.hal.science/tel-03219977v2

Submitted on 19 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Ensemble Algorithms and Analytic Combinatorics
in RNA Bioinformatics and Beyond

Yann Ponty

An Habilitation à Diriger des Recherches thesis of Paris-Saclay Université,
defended on May 22nd 2020 and evaluated by:

Cédric CHAUVE, Simon Fraser University, Canada

Léonid CHINDELEVITCH, Imperial College London, UK

Johanne COHEN, CNRS/Université Paris-Saclay

Philippe DUCHON, Université de Bordeaux

Ivo HOFACKER, University of Vienna, Autriche

Daniel MERKLE, University of South Denmark

Cyril NICAUD, Université Paris-Est, Marne la Vallée

Peter STADLER, University of Leipzig, Germany

## Abstract

Predictive Bioinformatics represents a major field of applications for combinatorial optimization techniques. Very often, an ensemble perspective, which not only consider the optimal solution but also fully embraces the set of suboptimal solutions, needs to be adopted. In this *Habilitation à Diriger des Recherches*, I present a series of algorithmic and combinatorial contributions, inspired by problems and questions arising in the study of RiboNucleic Acids (RNAs), in particular pertaining to their structural properties at the thermodynamic equilibrium.

I first describe a collection of generic and applied algorithmic techniques, enabling the efficient computation of statistical properties within search spaces. Such computations can be exact, or rely on unbiased estimates produced, using constrained sampling strategies, and are founded on a combinatorial (re-)interpretation of dynamic programming schemes. I then adopt a purely combinatorial point of view over search spaces, and establish asymptotic properties of classes of discrete objects arising in Bioinformatics, showcasing the unreasonable power of (a subset of) analytic combinatorics. Finally, I conclude with a collection of algorithmic results, obtained through the application of a wide array of techniques, in the context of RNA design, a field focused on combinatorial problems that are at the same time original, difficult, and relevant to the modern goals of biology and medicine.
# Contents

## 1 Introduction 5

1.1 Predictive Bioinformatics and its foundations 5

1.2 A crash course into RNA structure prediction 7

1.2.1 RNA structure 8

1.2.2 RNA folding prediction models and paradigms 9

1.2.3 The secondary structure 11

1.2.4 RNA 2D structure prediction from thermodynamic principles 13

1.3 Outline 16

## 2 Ensemble dynamic programming: techniques and analyses 17

2.1 A formal framework and basic algorithms 19

2.1.1 Dynamic programming as a rewriting system 20

2.1.2 Search space and suitability for ensemble applications 22

2.1.3 Classic optimization 24

2.2 Exact computation of Ensemble properties 25

2.2.1 Computing the partition function. 26

2.2.2 Probabilities in Boltzmann-Gibbs distributions (inside-outside) 27

2.2.3 Ensemble centroid and Maximum Expected Accuracy (MEA) 29

2.2.4 General moments of additive scores [150] 31

2.2.5 Classified DP with the Discrete Fourier Transform (DFT) [173, 174] 34

2.3 Probabilistic estimates 37

2.3.1 Foreword: On the number of samples [165] 38

2.3.2 Statistical sampling 39

2.3.3 Non redundant sampling [118, 133, 165] 41

2.3.4 Adaptive sampling of constrained sequences [18, 90, 158] 45

2.4 Applications of ensemble dynamic programming 48

2.4.1 (Recursive) simple type pseudoknots and kissing hairpins [150] 49

2.4.2 Dual partition functions and evolutionary robustness of RNAs 51

2.4.3 Unambiguous tree alignments [32, 33] 53

## 3 RNA design 56

3.1 Why do we design RNAs? 58

3.1.1 The potential of RNA design for statistical evolutionary studies 59
<table><tr><td>3.1.2</td><td>The different flavors of RNA design</td><td>59</td></tr><tr><td>3.2</td><td>Exact combinatorial negative design [85, 86]</td><td>61</td></tr><tr><td>3.2.1</td><td>Basic results</td><td>62</td></tr><tr><td>3.2.2</td><td>Design as a tree coloring problem</td><td>63</td></tr><tr><td>3.2.3</td><td>Structure-approximating design</td><td>65</td></tr><tr><td>3.3</td><td>Stochastic positive design under constraints</td><td>66</td></tr><tr><td>3.3.1</td><td>Sampling from the dual-partition function [117, 158, 160]</td><td>67</td></tr><tr><td>3.3.2</td><td>Avoiding and forcing motifs [210]</td><td>69</td></tr><tr><td>3.3.3</td><td>Multiple structures and extended design principles [88, 90]</td><td>71</td></tr></table>

4 Constrained random generation through rejection and beyond 73

<table><tr><td>4.1</td><td>It's not you, it's me: The unreasonable power of rejection</td><td>76</td></tr><tr><td>4.1.1</td><td>Correcting a (bounded) bias</td><td>77</td></tr><tr><td>4.1.2</td><td>Complexity aspects</td><td>78</td></tr></table>

<table><tr><td>4.2</td><td>Rejecting in dimension one and then some...</td><td>81</td></tr><tr><td>4.2.1</td><td>(Combinatorial) Boltzmann sampling</td><td>81</td></tr><tr><td>4.2.2</td><td>Multidimensional Boltzmann sampling</td><td>82</td></tr><tr><td>4.2.3</td><td>Beyond simplicity</td><td>88</td></tr></table>

<table><tr><td>4.3</td><td>Applications of multidimensional sampling</td><td>89</td></tr><tr><td>4.3.1</td><td>Controlling the GC-content [88, 90, 158, 160, 192, 193]</td><td>89</td></tr><tr><td>4.3.2</td><td>Dinucleotides content of protein-coding RNAs [209]</td><td>90</td></tr></table>

<table><tr><td>4.4</td><td>Redundancy in weighted sampling and countermeasures</td><td>92</td></tr><tr><td>4.4.1</td><td>Collisions in weighted languages [58, 78]</td><td>92</td></tr></table>

5 Asymptotic properties of RNA secondary structures and other trees 96

<table><tr><td>5.1</td><td>Basic tools in enumerative and analytic combinatorics</td><td>97</td></tr><tr><td>5.1.1</td><td>Context-free languages.</td><td>97</td></tr><tr><td>5.1.2</td><td>Enumeration and generating functions.</td><td>99</td></tr><tr><td>5.1.3</td><td>Basic singularity analysis</td><td>101</td></tr><tr><td>5.1.4</td><td>Useful extensions and shortcuts</td><td>105</td></tr></table>

<table><tr><td>5.2</td><td>Asymptotic combinatorics of RNA secondary structures</td><td>111</td></tr><tr><td>5.2.1</td><td>Expected 5'-3' distance [38]</td><td>113</td></tr><tr><td>5.2.2</td><td>RNA network properties [180]</td><td>115</td></tr><tr><td>5.2.3</td><td>RNA shapes: Abstract representations of structures [120]</td><td>119</td></tr><tr><td>5.2.4</td><td>Enumerating designable structures [206]</td><td>123</td></tr></table>

6 Conclusion and perspectives 127

Bibliography 129
## Foreword

The French academic system has enjoyed a long and eventful history, leading to many charming peculiarities. One such specificity is the requirement of junior tenured scientists to successfully defend an *Habilitation à Diriger des Recherches* (HDR) before being allowed supervise doctoral candidates in an official capacity. As part of the requirements for such a defense, the candidate produces a manuscript summarizing the research results obtained since his PhD. This is my contribution towards applying for an HDR in Computer Science at Université Paris-Sud. The individual nature of the evaluation dictates the use of the first-person singular, but an overwhelming majority of those results hinge critically on contributions from a community of collaborators, which I have attempted to acknowledge profusely.

The structure, ambition, scope and length of an HDR vary substantially between disciplines, institutions, and follow personal aesthetics. Some will elicit to contribute a general introduction for a collection of articles, while others will interlace past results with novel research material, trying to use the opportunity to outline general theories. Some will embrace the diversity of their past contributions, while others will strive to provide a unified perspective over research projects spanning one to two decades.

My ambition for this document is, by formalizing some of the intuitions underlying my past contributions, to reveal their triviality and, more seriously, contribute firm foundations for the systematic design and analysis of algorithms in RNA Bioinformatics and beyond. My main focus is therefore on the description of techniques, either classic in other fields or contributed over the course of my young career as a scientist, showcasing their power and level of generality by a brief description of their application in the context of RNA bioinformatics.

For these reasons, and partly because compact does not necessarily mean simple, I apologize in advance to the casual reader for explanations and digressions which may appear, at times, unnecessarily technical. It is unfortunately the price to pay to achieve my intended goal of using this document to lay out stable and explicit foundations for futures contributions to Bioinformatics.
# Chapter 1

## Introduction

Bioinformatics as a field of study is the poster child of **interdisciplinarity**. It is unified by an overarching objective to automate the processing and analysis of biological data, and plays an essential part in the production of knowledge in modern **Biology**. In the context of **Molecular Biology**, it is informed by models stemming from **Biophysics** and **Chemistry**, quantitatively parameterized by computational methods in **Statistics** and **Probability theory (Machine Learning)**, or analyzed at a theoretical level using **Theoretical Physics** and **Discrete Mathematics** techniques. Such quantitative models, typically in conjunction with parsimony arguments, represent the foundations of **predictive algorithmic methods**.

### 1.1 Predictive Bioinformatics and its foundations

**Predictive Bioinformatics** strives to produce methods which, from empirical data, predict phenomena far beyond our capacity for direct observations. Examples of such limitations include events that occur at the nanometer scale, in the distant past, or whose observation significantly interferes with (e.g. kill) a living system of interest. Given an established model, consistent with current biological knowledge, **predictive methods in bioinformatics** embrace some notion of search space induced by the input data. They elect one or several element(s) of the search space by maximizing some notion of score, usually analogous to a probability, provided by the model. The produced solution constitutes a *best bet* under (possibly implicit) assumptions of the model. Such methods are typically trained on reference data sets, for which a **ground-truth** is known, to calibrate a widely-varying set of parameters, and **validated** on independent data sets as an empirical test of their accuracy and, importantly, capacity of generalization.

Once a predictive method has been deemed satisfactory, a **leap of faith** occurs, following which the model and method are treated as uniformly correct. The subsequent predictions are considered as reflective of Nature itself, and treated as primary data in
<figure><img src="image_4.png" /><figcaption>Figure 1.1: **Rationale of ensemble analyses.** In combinatorial optimization, near-optimal solutions can be **representative** of the ensemble (A), with suboptimals concentrating around the optimal. Ensembles may also be fragmented (B), i.e. near-optimal solutions regroup into clusters whose representatives may differ from the optimal with respect to some features F<sub>A</sub> and F<sub>B</sub> of interest. The optimal may then be an **outlier**, and should be **treated with caution**, or even **disregarded** in favor of the ensemble centroid.</figcaption></figure>

the formulation of hypotheses and functional models in Biology. This abrupt change of perspective, treating a method with uttermost paranoia during the validation phase and with unabashed confidence in subsequent analyses, may appear dangerously optimistic at first sight. In fact, it is usually the root cause of the *cultural clash* experienced by scientists transitioning from exact sciences to Bioinformatics. However, Biology as field has evolved a **remarkable robustness** to misleading observations and appears, sometimes through impassioned controversies, to overcome biased observations caused by imperfect experimental devices, among which Bioinformatics methods can now be counted.

Indeed, predictions produced by Bioinformatics methods are typically falsifiable, i.e. they can be contradicted by future produced experimental data. This partly justifies the popularity of **exact algorithms** in molecular biology. Indeed, in the context of an exact optimization, failure to predict can be interpreted as a shortcoming of the underlying quantitative model, motivating and directing further refinements. By contrast, ill-defined optimization schemes (e.g. machine learning or heuristics), convolve possible modeling errors with optimization errors. Failure to predict does not necessarily point towards flaws of the model, and does not allow to refine them, as an improved model may induce worse optimization performances (and vice-versa, a broken clock being right twice a day), missing an opportunity for future developments.

Predictive optimization-based methods, however, suffer from a severe shortcoming in the
context of combinatorial search spaces. Indeed, they ultimately output a single solution, presumed to be optimally-likely in some sense. However, the dominating nature of the proposed solution may show limited robustness to, even modest, changes of the objective function, whose parameters are usually learned from data and thus subject to experimental errors. Moreover, an optimal solution may be poorly representative of the subspace of near-optimal solutions, so that it may be more advisable to consider a suboptimal – representative – solution rather than an optimal – non-representative – one.

These observations constitute the main motivations behind the development of **ensemble methods** in combinatorial Bioinformatics, treating the whole search space as an object of study rather than merely solving a *needle in a haystack* optimization problem. To that purpose, a probability distribution is postulated over the search space, and probabilistic analyses are performed. The optimal solution may remain relevant, but ensemble methods allow to go further, and assess notions of **support** for any given solution, or how distant features of the optimal solution are from the **centroid of the ensemble**. Indeed, as shown in Figure 1.1, ensembles associated with a given instance may either be concentrated around the optimal, or fragmented, leading to multiple clusters of diverse solutions, and multimodal distributions for features of interest.

## 1.2 A crash course into RNA structure prediction

The Bioinformatics of **RiboNucleic Acids** (RNAs) represent a very natural realm of application, and a source of constant inspiration, for ensemble analyses.

RNA constitutes a category of biomolecules, abstracted as **sequence of nucleotides** Adenine (A), Cytosine (C), Guanine (G) and Uracil (U), initially **transcribed** from a DNA template and further processed before reaching their cellular environment. They can form stable complexes with proteins, but also DNA and other RNAs, allowing them to **regulate genetic expression**. They can also perform **enzymatic functions**, i.e. process other molecules (or themselves), act as biosensors by undergoing conformational changes upon binding with small metabolites, and store the **entire genetic material** of certain viruses (e.g. HIV, SARS, 2019 nCoV).

Their dual capacity to store and process information is unmatched, leading current theories in evolution to consider RNA as the most likely candidate at the **origin of life** [93]. Such a versatility, illustrated in Figure 1.2, not only stems from the combinatorial nature of RNA sequences but also from its capacity to adopt one or several well-defined **structures**, driving the specificity of its interactions with other actors of the cellular world.
<figure><img src="image_5.png" /><figcaption>A – Basic RNA functions</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>B – Evolution of RNA families</figcaption></figure>

**Figure 1.2: Beyond coding for proteins, a growing list of RNA functions.** The versatility of RNA (A) is confirmed by the sustained growth (B) of the number of functional families in the RFAM database [84, 106], a majority of which are not coding for proteins, and associated with a conserved well-defined secondary structure.

### 1.2.1 RNA structure

Following Tinoco and Bustamante [184], RNA is believed to **fold** in a hierarchical fashion where **canonical base pairs**, mediated by hydrogen bonds, initially form a tree-like architecture called the **secondary structure**. Later in the folding process, the secondary structure is completed by complex topological motifs, including **pseudoknots** and **non-canonical base pairs** and motifs [114, 115]. Due to its combinatorial nature, the secondary structure is found at the core of successful computational methods for predicting the 3D structure of RNA [131]. It also supports the discovery of new functional families of **non-coding RNAs**, RNAs that do not (only) encode a protein, but (also) act directly on their environment through participation in catalytic and/or regulatory functions. A conserved **consensus secondary structure** is indeed a central object for many families in the RFAM database [84, 106], and a crucial element of the design of covariance models [62] used to find homologous RNAs across sequenced genomes.

At a 3D level, structure modeling is tackled using a mixture of comparative methods [134] and low-throughput/high resolution experimental techniques such as X-ray crystallography, Nuclear Magnetic Resonance (NMR), or Cryogenic Electron Microscopy (Cryo-EM). However, those techniques are extremely time-consuming, and require a preparation of molecules that may either not be feasible for certain RNAs, or can interfere with the folding process itself. These difficulties induce a **growing gap** between the number of functional families (±3000 as of Jan 2020), many of which requiring the adoption of a specific structure, and the number of families with experimentally-resolved 3D model available for at least one of their members (99 as of Jan 2020).

Recently, high throughput/low resolution alternatives have been proposed in the form
of improved **chemical probing** protocols, notoriously including SHAPE probing [47, 176, 198]. Those methods expose RNA to a chemical reagent, whose affinity towards individual nucleotides depend on the structure (and therefore partly reveal it, albeit in a stochastic and highly noisy fashion) and can be quantified using DNA and RNA sequencing technologies. The end-result of those methods are **1D reactivity** profiles that are not sufficient to fully characterize a structure, but greatly informative for further (computational) modeling.

### 1.2.2 RNA folding prediction models and paradigms

For all the aforementioned reasons, computational **structure prediction** methods are very relevant to the current objectives and challenges of RNA Bioinformatics. In the context of molecular biology, RNA structure prediction is mainly concerned with the prediction of one (or several) functional fold(s) for a given molecule. Since the influence of the cellular environment on RNA folding is difficult to fully characterize, and even harder to capture computationally, popular *in silico* methods adopt an approach inspired by **statistical mechanics**. They focus on stable conformations, having **low free-energies**, under the rationale that unstable structures are unlikely to be recognized by their partners, and play a reliable role in important phenotypic effects.

This notion can be formalized by considering the **Boltzmann-Gibbs distribution**, where any possible structure S for a given RNA w is observed with probability

$$P_w(S) \propto e^{-E(w,S)/RT}, \qquad (1.1)$$

where $E(w, S)$ represents the free-energy in kcal.mol⁻<fnref n="1" /> of the (w, S) pair, R the Boltzmann constant (1.9872 10⁻³ kcal.mol⁻<fnref n="1" />.K⁻<fnref n="1" />) and T the absolute temperature in Kelvin. This distribution can be thought as the **stationary distribution** of a **continuous time Markov chain**, illustrated in Figure 1.3 on a toy example. In this process, a transition between two conformations is chosen with probability/rate that only depend on the energy difference, thus respecting a **detailed balance**. Starting from the fully-unpaired/empty structure<fnref n="1" />, having initial probability 1, the probability mass diffuses within the space of conformations, and ultimately reaches the **thermodynamic equilibrium** when probabilities no longer evolve with time. The evolution with time of probabilities/concentrations can be determined using efficient numerical integration [201], as illustrated in Figure 1.3.B.

This vision, inspired by statistical mechanics, is at the core of **historical paradigms** in the field of RNA structure prediction:

<footnote n="1">¹More sophisticated models may consider a co-transcriptional folding [111] of the nascent transcript, leading to a non-binary initial distribution at the beginning of the process.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 1.3: RNA folding paradigms. The process of RNA folding can be abstracted as a Continuous Time Markov Chain (CTMC – A), over a discrete state space consisting of (a subset of) the secondary structures adopted by an RNA. Starting from an initial distribution, typically assigning full probability to the open chain, kinetic studies attempt to characterize the evolution of the probability distribution with time (B). At the thermodynamic equilibrium, the stationary distribution of the CMTC is reached and probabilities cease to evolve with time. The most probable structure is then the most stable one, i.e. the MFE structure. However, RNA degradation may occur before the equilibrium is reached, and the dominant structures at finite time may represent more promising candidates for functional hypotheses.</figcaption></figure>

**MFE paradigm.** Early computational methods for structure prediction [95, 143, 213] focused their effort on producing (one of) the optimally stable, or **Minimum Free-Energy (MFE)**, structure(s), under the rationale that the MFE structure has highest probability at the thermodynamic equilibrium;

**Boltzmann ensemble paradigm.** While optimal, the probability of the MFE structure alone can be abysmally small, and decreases exponentially with the RNA length. The MFE structure can therefore be **overwhelmed**, at the thermodynamic equilibrium, by a set of alternative structures having comparable stability and, possibly, very different characteristics. This motivates studies of the **expected properties** in the Boltzmann-Gibbs distribution, hinging on an efficient computation of the **partition function** [130]

$$Z = \sum_{S} e^{-E(w,S)/RT}.$$

When these properties diverge from those of the MFE, it may be more relevant to consider structures that are more accurate representatives of the ensemble of structures. One such structure is the **Maximum Expected Accuracy** structure [55, 121], the structure minimizing some notion of expected distance to the rest of the ensemble. Alternatively a statistically representative, Boltzmann distributed, set of structures can be randomly generated, clustered at a structural level to identify
alternative conformations, and eliminate outliers. A **centroid structure** can then be elected for each cluster [52];

**The kinetics paradigm.** More recently, there has been a growing awareness of the importance of **out-of-equilibrium effects** in the function(s) carried out by RNA. Indeed, in a cellular context, RNA is constantly transcribed and degraded by enzymes. Depending on the precise dynamics of these concurrent processes, a population of RNAs may simply be degraded before reaching its stationary distribution, so the thermodynamic paradigm may not provide an accurate picture of the functional structure. This generally applies to RNAs whose **energy landscapes** feature substantial **barriers**, leading to a slow convergence towards the equilibrium. Kinetics analyses are thus concerned with the structural behavior of RNA *before* reaching the equilibrium.

Examples of kinetics effects are suspected to include **riboswitches**, bistables RNAs which are observerd in *on* and *off* conformations in the presence/absence of a **ligand**, a small molecule, believed to have insufficient contribution to the free-energy to invert the relative stabilities of the *on* and *off* states. Current models are thus based on kinetics, and postulate that the ligand modifies an energy barrier, leading to a faster/slower convergence towards the thermodynamic equilibrium [65]. More generally, **co-transcriptional folding**, the folding of RNA during transcription, reveals the importance of kinetic effects. Indeed, this phenomenon would be without effect at the thermodynamic equilibrium, since the stationary distribution of a (ergodic) Markov chain does not depend on its initial distribution, so it would should not matter at the equilibrium.

In the rest of this document, I will mainly focus on algorithmic strategies relevant to the MFE and Boltzmann equilibrium paradigms, although certain were have designed with kinetics analysis in mind [133, 180]. Kinetics analyses are indeed much more time-consuming, and are associated with computational problems that are routinely NP-hard [126]. Efficient heuristics and methods for analyzing kinetics are, however, the object of ongoing projects within the RNA Bioinformatics community.

### 1.2.3 The secondary structure

A **RNA secondary structure** of size n represents the outcome of a folding process, and focuses on a subset of **base-pairs**, mediated by hydrogen bonds. For essentially computational reasons [3, 123, 175] this definition forbids crossing base pairs, also called **pseudoknots** due to their ability to induce complex topologies [205]. Moreover, any nucleotide can only be involved in a single base pair, since additional partners would involve non-canonical edges [114]. Finally, most definitions rule out base pairs between
<figure><img src="image_8.png" /><figcaption>Figure 1.4: Various representations for RNA secondary structures.</figcaption></figure>

proximate positions, due to **steric effects** inducing geometric constraints, leading to a **minimal number θ** of unpaired positions between paired positions.

**Definition 1.2.1 (RNA secondary structure):** An RNA secondary structure S of length n is a set of base-pairs (i, j), 1 ≤ i < j ≤ n, such that:

* Each position is **monogamous**, ∀(i, j) ≠ (i', j') ∈ S : {i, j} ∩ {i', j'} = ∅;

* **Minimal distance θ** between paired nucleotides, ∀(i, j) ∈ S : j - i > θ;

* **No pseudoknot** allowed, ∀(i, j), (i', j') ∈ S, i < i' : (j' < j) or (j < i').

The **conformational space** associated with a sequence of length n is simply $S_n$, the set of secondary structures over n nucleotides.

In typical *ab initio* RNA structure prediction problems, the input is a sequence of nucleotides $w \in \{A, C, U, G\}^*$. This space of secondary structure is then usually restricted to structures consisting of **canonical base pairs** $B := \{{\{G, C\}, \{A, U\}, \{G, U\}\}}$. In other words, for any **valid secondary structure** $S \in S$ one has:

$$ \forall (i, j) \in S : \{w_i, w_j\} \in B. $$

The secondary structure can be drawn in a variety of ways, as illustrated by Figure 1.4, many of which being supported by our popular software VARNA, developed in collaboration with Kevin Darty and Alain Denise [44].

The stability of a secondary structure is assessed using a free-energy model. The most popular such model is the **Turner nearest-neighbor free-energy model** [185], which associates experimentally-determined **free-energy** contributions to structural motifs, called **loops** (see Figure 1.5). The energy of any given secondary structure is then additively defined, i.e. obtained by summing the contributions of the various loops appearing in the structure.
<figure><img src="image_9.png" /><figcaption>Figure 1.5: Loop decomposition supporting the Turner nearest neighbor model [185]. Free-energies are associated to each loop type, precise topology and content in nucleotides, determined by, or extrapolations from, direct experimental measurements.</figcaption></figure>

### 1.2.4 RNA 2D structure prediction from thermodynamic principles

The secondary structure provides a convenient discretization of the conformational space, and allows a reformulation of several folding paradigms into combinatorial problems. For instance, predicting the most stable secondary structure (MFE paradigm) amounts to finding the secondary structure having minimal free-energy according to a chosen model. Since the overall and expected number of secondary structure of length n both grow exponentially on the length [194, 212], a proof being provided in Section 5.2, brute-force optimization is not realistically feasible.

However, as illustrated in Example 1 below, the tree-like nature of secondary structure, and the independence of contributions within the energy model, make those optimizations amenable to a $\Theta(n^3)$ **dynamic programming** (DP), as initially shown by Nussinov and Jacobson [143] in a simplified energy model. The algorithm was later extended to capture loops and additional features of the Turner free-energy model by Zuker and Stiegler [213], leading to a $\Theta(n^4)$ time algorithm (empirically running in $\Theta(n^3)$ time). This algorithm was then extended by [204] to generate the exhaustive collection of secondary structures within a free-energy range $\Delta$ of the MFE structure.

#### Example 1: RNA BP folding – Running example

A historical problem, first efficiently solved by Nussinov and Jacobson [143] is the prediction of the **MFE structure in a simple base pairing model**, i.e. having maximum number of pairs, with only nucleotides in B being allowed to pair. Our goal is thus to compute some $S^* \in S_n$ such that

$$|S^*| = \max_{S \in S_n \text{ such that } V(i,j) \in S, \{w_i, w_j\} \in B} |S|.$$

A classic **dynamic programming scheme**, inspired by Waterman [194], considers each interval $[i, j] \subset [1, n]$ and their associated optimal structure $S^*_{i,j}$ for $[i, j]$:
$$ \text{\textasciitilde{i}} \text{\textasciitilde{j}} = \text{\textasciitilde{i}} \text{\textasciitilde{i+1}} \text{\textasciitilde{j}} + \text{\textasciitilde{i}} \text{\textasciitilde{i+1}} \text{\textasciitilde{k-1}} \text{\textasciitilde{k}} \text{\textasciitilde{k+1}} \text{\textasciitilde{j}} $$

Reasoning on the fate of the first position i within $S_{i,j}^*$, one only needs to consider three cases:

1. $i + \theta \ge j$, so no base pair can be formed, and we have $S_{i,j}^* := \emptyset$;

2. $i + \theta < j$ and i is left unpaired, so any base pair in $S_{i,j}^*$ involves positions in $[i+1, j]$, and we have $S_{i,j}^* := S_{i+1,j}^*$;

3. $i + \theta < j$ and i is paired to some $k \in [i+1, j]$ in $S_{i,j}^*$ subject to $\{w_i, w_k\} \in \mathcal{B}$. Beyond $(i, k)$, any base pairs in $S_{i,j}^*$ involves positions in $[i+1, k-1]$ in $[k+1, j]$ (but not both since pseudoknots are forbidden), so we get $S_{i,j}^* := \{(i, k)\} \cup S_{i+1,k-1}^* \cup S_{k+1,j}^*$.

This immediately suggests a **recurrence for the maximum number of base pairs** $E_{i,j}$, achieved by some secondary structure for the region $[i, j] \in [1, n]$ of a given RNA w:

$$ \begin{array}{l} i + \theta \ge j : E_{i,j} := 0, \\ i + \theta < j : E_{i,j} := \max \left\{ \begin{array}{ll} E_{i+1,j} & (\text{i unpaired}) \\ \max_{k=i+\theta+1}^{j} 1 + E_{i+1,k-1} + E_{k+1,j} & (\text{i paired}) \end{array} \right. \end{array} $$

Note that the computation of $E_{i,j}$ only requires values $E_{i',j'}$ such that $\| [i', j'] \| < \| [i, j] \|$. A **dynamic programming algorithm** for computing $|S^*| := E_{1,n}$ will consider intervals $[i, j] \subseteq [1, n]$, **ordered by increasing span** $[[i, j]]$, and compute $E_{i,j}$ using precomputed $E_{i',j'}$, storing the result before proceeding to a (possibly) larger interval.

The **correctness** of the algorithm hinges on the fact that, while computing over $E_{i,j}$, all $[i', j']$ such that $\| [i', j'] \| < \| [i, j] \|$ are already been processed, so that $E_{i',j'}$ in the right hand side can be used. The number of different $E_{i,j}$ terms grows in $\Theta(n^2)$. Moreover, each $E_{i,j}$ can be computed in $O(n)$ arithmetic operations, assuming that all $E_{i',j'}$ in the right-hand-side are accessible in $O(1)$ operations. The **overall complexities** of the algorithm are therefore in $\Theta(n^3)/\Theta(n^2)$ time and space.

Once the maximum number of base pairs is computed, one can **backtrack** through these recursions, starting from $[1, n]$, by figuring out at each step (one of) the case(s) that contributed to the max. Concretely, one will define a recursive function Backtrack such that Backtrack(i, j) := ReturnØ when $i + \theta \ge j$ and, for $i + \theta < j$:

$$ \text{Backtrack}(i, j) := \begin{cases} \text{if } E_{i,j} = E_{i+1,j} : & \text{return Backtrack(i,j)} \\ \text{if } E_{i,j} = 1 + E_{i+1,k-1} + E_{k+1,j}, k \in [i + \theta + 1, j] \text{ and } \{w_i, w_k\} \in \mathcal{B} : \\ & \text{return}\{(i, k)\} \cup \text{Backtrack}(i, k-1) \cup \text{Backtrack}(k+1, j) \end{cases} $$

Note that the two cases are not necessarily mutually-exclusive, i.e. multiple structures may achieve maximal number of base pairs, and an arbitrary-chosen optimal one will be returned.

A second breakthrough came with the advent of **comparative folding**, pioneered by the work of Sankoff [167], who automated the practices of manual modelers [134], and proposed an algorithm for the **simultaneous folding and alignment** of RNA. While associated with extreme computational demands ($\Theta(n^6)$ for aligning two homologous sequences), its impressive gain in predictive power and reliability compared to *ab initio* techniques led it to be adapted within many subsequent methods, undergoing active
developments [128, 179, 199, 200].

A revolution came with the McCaskill [130] algorithm, which championed a transition towards the **Boltzmann ensemble** paradigm. Namely, McCaskill observed that the **partition function**

$$Z = \sum_{S} e^{-E(w,S)/RT}$$

could be computed in essentially $\Theta(n^3)$ time, through a simple change of algebra to the dynamic programming equations of [213]. Moreover, he showed that the same decomposition could be adapted into an instance of the inside/outside algorithm[10, 113], leading to a $\Theta(n^4)$ algorithm for computing the **base-pair probabilities** within the Boltzmann ensemble. Such **ensemble properties** provide a notion of support for structures predicted in the MFE paradigm [127], but also allow to predict a structure that more adequately represents the ensemble than the MFE, such as the **Maximal Expected Accuracy** structure [55, 121].

The flexibility, and scope of applications, of ensemble methods was greatly extended by the contribution by Ding and Lawrence [51] of a **stochastic backtrack** algorithm, allowing to produce a random, Boltzmann distributed, secondary structure in time $\Theta(n^2)$. This enabled the implementation of **statistical estimators** for many features, including some that cannot be captured by dynamic programming. This also paved the way for a statistical estimation of the **dominant conformations** within the Boltzmann ensemble, using a combination of sampling and clustering [52], used by subsequent methods [177].

Finally, a strong emphasis has been recently put on the development of **integrative methods** that exploit the availability of (partial/noisy) experimental information. In the early 2000s, enzymatic and chemical probing data were integrated as **hard constraints** by Mathews et al. [129] within DP equations for RNA structure prediction. Following the development of quantitative experimental methods, such as the SHAPE technologies [176, 198], probing data are now incorporated as **soft constraints**, a.k.a. **pseudo-energies** within DP schemes for prediction. Such energy terms can be thought as *shifting* the Boltzmann ensemble towards areas of greater compatibility with the **reactivity profiles** produced by the probing experiments.
## 1.3 Outline

In this manuscript, I attempt to provide an unifying view on the algorithmic and theoretical concepts, used in a series of personal (frequently collaborative) contributions in Bioinformatics, Theoretical Computer Science and Discrete Mathematics. Those contributions are of different nature, and address a variety of objects of study, yet they share a consideration of ensembles of combinatorial objects, and are based on an interpretation of enumerations schemes as algorithmic principles.

In Chapter 2, I reformulate **ensemble dynamic programming** algorithms, contributed in the context of Bioinformatics within a unified framework. My contributions draw on a strong connection between (ensemble) dynamic programming and enumerative combinatorics, and put a strong emphasis on the design of dynamic programming schemes, whose productions are bijectively associated with an ensemble of objects of interest.

Chapter 3 focuses more specifically on algorithmic methods for RNA design, the algorithmic construction of new RNA sequences achieving a certain function, here through the adoption of a given structure. In this context, I revisit the classic inverse folding of RNA at under (essentially) a base pair maximization model, and show families for design can be approximated in an original sense. I also develop a vision inspired by random generation for the positive design problem, where one attempts to favor the affinity towards a given fold.

Chapter 4 describes rejection-based techniques and algorithms for a controlled random generation of combinatorial objects. Again, the initial focus is to present those techniques in an application-agnostic setting, later to describe some of their applications in Bioinformatics. Finally, I mention some analyses of the, arguably uninformative, redundancy within samples and the shortcomings of a rejection-based strategy to overcome it.

Finally, Chapter 5 summarizes a series of analyses focusing on the asymptotic properties of objects occurring in (RNA) bioinformatics. Those include RNA secondary structures, for which a careful application of analytic combinatorics principles allow to derive properties in the homopolymer model, but also other objects with connection to the design and analysis of algorithms.
Chapter 2

Ensemble dynamic programming:
techniques and analyses

Following Bellman [14], a **dynamic programming algorithm** is usually stated as a system of recursive equations, relating the optimal value of an **objective function** over a certain argument, or problem, to a set of its values over other arguments, or sub-problems. Additionally, such an equation has induce an **acyclic computation**, meaning that the set of arguments that are used by following the recursive calls, can be totally ordered in a way that is consistent with the left-to-right transitions in the system. Such a definition overlooks an important aspect of a dynamic programming scheme: the semantics associated with the choice of one of the available left-to-right transitions.

In this manuscript, we adopt a perspective dually inspired by **enumerative combinatorics**, where the choice of a derivation provides partial information regarding the final solution. Such an enumerative perspective allows to define notions of unambiguity, correctness and completeness with respect to a given search space which, if fulfilled by a DP scheme, unlock a variety of algorithms to make statements regarding the **ensemble of solutions**. Such analyses go **beyond the optimization of an objective functions**, and attempt to answer selected questions involving the ensemble of candidate solutions:

* What is the support of an optimal solution?

• What is the number of near-optimal solutions?

* Are all near-optimal solutions similar? how diverse are they?

* What is the centroid solution, i.e. most similar to other near-optimal solutions?

* What are the average properties of near optimal solutions?

* Which distributions of properties are expected within near optimal solutions?

Such questions require the definition of a **probability distribution**, assigning proba-
bilities to elements of the search space increasingly with their value for the objective
function.
**Ensemble analyses**, whose first instances can probably be traced back to the origins of Natural Language Processing [10, 188], enjoy a great popularity in Bioinformatics. In particular, the early age of RNA bioinformatics saw seminal contributions from scientists with a strong combinatorial culture in Discrete Mathematics, Theoretical Physics and Computer Science, including Michael Waterman [194], Michael Zuker [212, 213], David Sankoff [167], Walter Fontana, Peter Stadler, Peter Schuster and Ivo Hofacker [74, 95, 96]....

The popularity of ensemble analyses in RNA Bioinformatics also stems directly from their origin in **statistical mechanics**, inspiring the seminal contribution of John McCaskill [130] who turned a popular DP scheme for energy minimization into an algorithm to compute the **partition function** of the Boltzmann ensemble. This central quantity allowed to derive expected properties of RNA at the **thermodynamic equilibrium**, starting with the **base pairing probabilities**, which have become central to modern methods for comparative folding of RNA [199, 200], based on a simultaneous folding and alignment of RNAs pioneered by David Sankoff [167].

In this context, my efforts have been focusing, on: i) the design of novel DP schemes amenable to ensemble analyses, associated with challenges of variable technicality; and ii) the development of new techniques to enable, or accelerate, ensemble analyses.

**Outline.** Section 2.1 describes a unifying framework for **ensemble dynamic programming**, used in Section 2.2 to reformulate classic and novel tools to perform an **exact computation of ensemble properties**. Section 2.3 focuses on **statistical sampling**, providing estimates for, possibly complex, Ensemble properties. Finally, Section 2.4 presents a selection of **algorithms and results** obtained in various areas of Bioinformatics using this general framework.

The following summarizes, and attempts to unify, contributions described within the following list of articles published in journals and/or presented at conferences.

**Associated contributions**

J. Waldispühl and Y. Ponty. An unbiased adaptive sampling algorithm for the exploration of RNA mutational landscapes under evolutionary pressure. **Journal of Computational Biology**, 18(11):1465–79, Nov. 2011

<figure><img src="image_10.png" /></figure>

J. Waldispühl and Y. Ponty. An unbiased adaptive sampling algorithm for the exploration of RNA mutational landscapes under evolutionary pressure. In **RECOMB 2011**, volume 6577 of *Lecture Notes in Computer Science*, pages 501–515, Vancouver, Canada, Mar. 2011. Springer Berlin / Heidelberg

Y. Ponty and C. Saule. A Combinatorial Framework for Designing (Pseudoknotted) RNA Algorithms. In **WABI 2011**, Saarbrucken, Germany, 2011

S. Sheikh, R. Backofen, and Y. Ponty. Impact Of The Energy Model On The Complexity Of RNA Folding With Pseudoknots. In **CPM 2012**, volume 7354 of *Combinatorial Pattern Matching*, pages 321–333, Helsinki, Finland, July 2012. Juha Kärkkäinen, Springer
P. Rinaudo, Y. Ponty, D. Barth, and A. Denise. Tree decomposition and parameterized algorithms for RNA structure-sequence alignment including tertiary interactions and pseudoknots. In WABI 2012, tba, Ljubljana, Slovenia, Sept. 2012. University of Ljubljana

E. Senter, S. Sheikh, I. Dotu, Y. Ponty, and P. Clote. Using the Fast Fourier Transform to Accelerate the Computational Search for RNA Conformational Switches. PLoS ONE, 7(12):e50506, Dec. 2012

<figure><img src="image_11.png" /></figure>

E. Senter, S. Sheikh, I. Dotu, Y. Ponty, and P. Clote. Using the Fast Fourier Transform to accelerate the computational search for RNA conformational switches (extended abstract). In RECOMB 2013, Beijing, China, Apr. 2013

V. Reinharz, Y. Ponty, and J. Waldispühl. Using Structural and Evolutionary Information to Detect and Correct Pyrosequencing Errors in Noncoding RNAs. Journal of Computational Biology, 20(11):905–19, Nov. 2013

<figure><img src="image_12.png" /></figure>

V. Reinharz, Y. Ponty, and J. Waldispühl. A linear inside-outside algorithm for correcting sequencing errors in structured RNA sequences. In RECOMB 2013, Beijing, China, Apr. 2013

C. Chauve, Y. Ponty, and J. P. P. Zanetti. Evolution of genes neighborhood within reconciled phylogenies: an ensemble approach. BMC Bioinformatics, 16(Suppl 19):S6, Dec. 2015

<figure><img src="image_13.png" /></figure>

C. Chauve, Y. Ponty, and J. P. P. Zanetti. Evolution of genes neighborhood within reconciled phylogenies: an ensemble approach. In BSB 2014, volume 8826 of Advances in Bioinformatics and Computational Biology, pages 49–56, Belo Horizonte, Brazil, Oct. 2014. Springer

A. Rajaraman, C. Chauve, and Y. Ponty. Assessing the robustness of parsimonious predictions for gene neighborhoods from reconciled phylogenies. In ISBRA 2015, volume 9096, pages 260–271, Norfolk, Virginia, United States, June 2015

E. Jacox, C. Chauve, G. J. Szöllösi, Y. Ponty, and C. Scornavacca. ecceTERA: Comprehensive gene tree-species tree reconciliation using parsimony. Bioinformatics, 32(13):2056–2058, July 2016

V. Reinharz, Y. Ponty, and J. Waldispühl. Combining structure probing data on RNA mutants with evolutionary information reveals RNA-binding interfaces. Nucleic Acids Research, 44(11):e104–e104, 2016

W. Duchemin, Y. Anselmetti, M. Patterson, Y. Ponty, S. Bérard, C. Chauve, C. Scornavacca, V. Daubin, and E. Trier. DeCoSTAR: Reconstructing the ancestral organization of genes or genomes using reconciled phylogenies. Genome Biology and Evolution, 9(5):1312–1319, 2017

J. Deforges, S. De Breyne, M. Ameur, N. Ulryck, N. Chamond, A. Saaidi, Y. Ponty, T. Ohlmann, and B. Sargueil. Two ribosome recruitment sites direct multiple translation events within HIV1 Gag open reading frame. Nucleic Acids Research, 45(12):7382–7400, July 2017

C. Chauve, J. Courtiel, and Y. Ponty. Counting, generating, analyzing and sampling tree alignments. International Journal of Foundations of Computer Science, 29(5):741–767, 2018

<figure><img src="image_14.png" /></figure>

C. Chauve, J. Courtiel, and Y. Ponty. Counting, generating and sampling tree alignments. In ALCOB 2016, volume 9702, pages 53–64, Trujillo, Spain, 2016. Springer

S. Hammer, W. Wang, S. Will, and Y. Ponty. Fixed-parameter tractable sampling for RNA design with multiple target structures. BMC Bioinformatics, 20(1):209, Dec. 2019

<figure><img src="image_15.png" /></figure>

S. Hammer, Y. Ponty, W. Wang, and S. Will. Fixed-Parameter Tractable Sampling for RNA Design with Multiple Target Structures. In RECOMB 2018, Paris, France, 2018

## 2.1 A formal framework and basic algorithms

In order to refactor both classic and novel dynamic programming-based algorithms, we introduce a formal framework. Largely inspired from Juraj Michalik's PhD [132], it can be seen as an operational version of the declarative proposal of Giegerich and Touzet [81], modeling dynamic-programming processes as inverse coupled term-rewriting systems [108]. Some of its features were previously introduced in collaboration with Cédric Saule [150], based on an oriented hypergraph framework pioneered by Finkelstein and Roytberg [66], independently pursued by Huang and Chiang in the context of natural language process-
ing [98]. It is also worth mentioning a substantial overlap with the book of Miklos [135] dedicated to the complexity of counting and sampling, of which our framework covers some of the easier (polynomially-solvable) cases.

## 2.1.1 Dynamic programming as a rewriting system

Denote by $w$ an **instance** (e.g. sequence, tree, graph...) of a combinatorial problem. An instance implicitly defines a discrete **universe** $U_w$, of which the **search space** $\Omega_w \subset U_w$ of a combinatorial algorithm is a subset. Next, we need to describe how elements of the search space are decomposed, or conversely generated, by recursive constructors.

**Definition 2.1.1 (Constructor):** A **constructor** of **arity** $k$ is a function $\lambda : U^k \to U$ that returns/creates a novel element the search space $\Omega$ from $k$ elements of $\Omega$.

In other words, a constructor is a function that **assembles** a candidate solution from a collection of (smaller) candidate solutions to subproblems. Constructors with zero **arity**, i.e. constant functions, are called **atoms** and constitute base cases in the classic recursive exposition of dynamic programming. To mark this distinction between a constructor $\lambda$ seen as a function, e.g. used to label the nodes of terms (see definition below), and its evaluation, we use the notation $\lambda \rightsquigarrow v$. Denote by $\Lambda$ the **set of all constructors** for a given DP scheme.

Prior to any definition of a dynamic programming scheme, one needs to introduce a **state space** $Q$. Any state $q \in Q$ represents a (sub)problem encountered while solving the recursive DP computation. The spirit of dynamic programming is to solve a given problem by solving a number of (smaller) problems, depending on the associated state and the instance. This dependency is materialized in our formalism by **derivations**, each decorated by a combinatorial **constructor**.

**Definition 2.1.2 (Derivation):** A **derivation** is a tuple

$$ (q, (q, \dots, q_k), \lambda) \in Q \times Q^* \times \Lambda, $$

denoted as $q \xrightarrow{\lambda} q_1, \dots, q_k$ such that:

* $q$ is the origin;

* $(q_1, \dots, q_k) \in Q^*$ is the production, i.e. an ordered list of states (a.k.a. subproblems) that have to be solved in order to solve $q$;

* $\lambda \in \Lambda$ is a constructor of **arity** $k$.
We can now define a (combinatorial) dynamic programming scheme as a system of equations, coupled with a derivation system.

**Definition 2.1.3 (Dynamic Programming Scheme):** A dynamic programming scheme $\Delta$ is a tuple $(Q, q_w, \delta)$, where:

* $Q$ is the state space;

* $q_w \in Q$ is the initial state;

* $\delta \subset Q \times Q^* \times \Lambda$ is an **acyclic**, i.e. non transitively self-referential, set of derivations.

The **acyclicity condition** forbids any state to (transitively) derive into itself, and is essential for algorithmic considerations. Note that, while the instance does not explicitly appear in this definition, its precise content is at the origin of the state space and lists of derivations. To illustrate this bundle of abstract definitions, let us reformulate our running example.

### Example 2: RNA BP folding – Formalized dynamic programming

The DP scheme $\Delta_{BP}$ for RNA 2D structure prediction problem, introduced in Section 1.2.3 can be formalized as follows: the **instance** is $w \in \{A, C, U, G\}^*$, a sequence of nucleotides; the **universe** is $U := \cup_{n' \le n} S_{n'}$; the **states** are the intervals $Q := \{[i, j] | 1 \le i \le j \le n\}$, $n := |w|$ of the input sequence $w$; and three types of **constructors** are sufficient to generate all structures:

* $\lambda_{\emptyset}^{[i,j]} \rightsquigarrow \emptyset \quad \rightarrow \quad \text{Atom, returning the empty secondary structure;}$

* $\lambda_{\bullet}^{[i,j]}(S) \rightsquigarrow S \quad \rightarrow \quad \text{Leaves position i unpaired;}$

* $\lambda_k^{[i,j]}(S, S') \rightsquigarrow S \cup S' \cup \{\{i, k\}\} \quad \rightarrow \quad \text{Adds base pair } \{i, k\} \text{ to two substructures } S \text{ and } S'$.

The **derivations** in $\delta$ consist of:

* Terminal derivations:

$$ \forall[i, j] \subset [1, n], j \le i + \theta \rightarrow (i, j) \xrightarrow{\lambda_{\emptyset}^{[i, j]}} \epsilon $$

Semantics: Position i left unpaired, no further processing required (sequence too short to support base pair), empty structure $\lambda_{\epsilon}$ returned;

* Unpaired derivations:

$$ \forall[i, j] \subset [1, n], j \le i + \theta : (i, j) \xrightarrow{\lambda_{\bullet}^{[i, j]}} (i + 1, j) $$

Semantics: Position i is left unpaired, requires processing of interval $[i+1, j]$, (optimal) structure built over $[i+1, j]$;

* Paired derivations:

$$ \forall[i, j] \subset [1, n], \forall k \in [i + \theta + 2, j], \{w_i, w_k\} \in \mathcal{B} : (i, j) \xrightarrow{\lambda_k^{[i, j]}} ((i + 1, k - 1), (k + 1, j)) $$

Semantics: Position i paired with k, requires processing of intervals $[i+1, k-1]$ and $[k+1, j]$, (optimal) structure built from $[i+1, k-1]$ and $[k+1, j]$, augmented with base pair $\{i, k\}$;

Derivation rules, in conjunction with constructors, define the search space explored (or, conversely, generated) by a DP scheme. To reason on the relationship between
constructors and elements of the search space, we make a subtle distinction between a **term** T, a tree-like hierarchy of constructors produced by a complete series of derivations, and its **evaluation** $\phi(T)$ as an element of the search space.

For instance, using the constructors and semantics of Example 2, we have

<figure><img src="image_16.png" /></figure>

## 2.1.2 Search space and suitability for ensemble applications

The terms space and search space produced by a given dynamic programming scheme can then be defined recursively as follows.

**Definition 2.1.4 (Terms of a DP scheme):** The **terms set** $T_{\Delta}$ generated by a dynamic programming scheme $\Delta = (Q, q_w, \delta)$ is defined as $T_{\Delta} := T_{q_w}$ where, for any $q \in Q$, one has:

$$ T_q = \bigcup_{\substack{\lambda \\ q \to q_1 \cdots q_k \in \delta}} \left\{ \bigwedge_{t_1 \ t_2 \cdots t_k}^{\lambda} \mid t_1 \in T_{q_1}, t_2 \in T_{q_2}, \dots, t_k \in T_{q_k} \right\}. $$

Terms represent the syntactical structure of elements of the search space, which we now define.

**Definition 2.1.5 (Search space of a DP scheme):** The **search space** $\Omega_{\Delta} := \Omega_{q_w}$ generated by a dynamic programming scheme can be similarly defined through

$$ \begin{aligned} \Omega_q &= \{\phi(T)\}_{T \in T_q} \\ &= \bigcup_{q \to q_1 \cdots q_k \in \delta} \{\lambda(s_1, \dots, s_k) \mid (s_1, \dots, s_k) \in \Omega_{q_1} \times \cdots \times \Omega_{q_k}\}. \end{aligned} $$

Equipped with these notions, we can now define the properties of a dynamic programming scheme, that will connect it to an underlying reality.
**Definition 2.1.6 (Completeness and unambiguity of a DP scheme):** A DP scheme $\Delta = (Q, q_w, \delta)$ is:

1. **Unambiguous** if and only if every element of the search space can be generated in only one way, i.e. $\phi$ is bijective between $\mathcal{T}_\Delta$ and $\Omega_\Delta$;

2. **Complete** with respect to a targeted search space $\Omega^*$ if and only if every element in $\Omega^*$ is considered by $\Delta$, i.e. one has $\Omega_\Delta = \Omega^*$.

These two notions are crucial for ensemble applications of application. Indeed, in combination, they allow to use a given DP scheme to extract relevant properties of a preexisting search space.

Let us now distinguish elements within the search spaces, by introducing **scoring functions** that will map numerical values with each constructors and, in turn, to terms/elements of the search space

**Definition 2.1.7 (Additive scoring function):** An **additive scoring function** $f: \Lambda \to \mathbb{R}$ associates a numerical value to each constructor, such that the **score** $f(T)$ of a term $T \in \mathcal{T}$ is defined as

$$f(T) := f(\phi(T)) := \sum_{\lambda \in \mathcal{T}} f(\lambda).$$

An important property of a dynamic programming scheme lies in its ability to **emulate a given function** defined over its search space, by using a suitable scoring of its constructors/derivations. Such function could represent an **objective function** in the context of an optimization, or help induce a desired probability distribution over the search space.

**Definition 2.1.8 (Correctness of a DP scheme):** Let $\Delta$ be a DP scheme, coupled with a scoring function $f: \Lambda \to \mathbb{R}$. A pair $(\Delta, f)$ is **correct**, with respect to a given function $F: \Omega_\Delta \to \mathbb{R}$, if and only if $F(\phi(T)) = f(T), \forall T \in \mathcal{T}$.

By extension, we say that a DP scheme $\Delta$ is **correct** with respect to a function $F$ if and only if there exists a scoring function $f$ such that $(\Delta, f)$ is correct.

### Example 3: RNA BP folding – Unambiguity/completeness/correctness

The **unambiguity** of $\Delta_{BP}$ can be proven by considering two terms $T$ and $T'$, $T \neq T'$. Consider the first position from the root where constructors $\lambda$ and $\lambda'$, such that $\lambda \neq \lambda'$ are found in $T$ and $T'$ respectively. Since their paths to the root encounter the same constructors, $\lambda$ and $\lambda'$ are of the form $\lambda_k^{[i,j]}$ or $\lambda_k^{[i,j]}$ (but not $\lambda_\emptyset^{[i,j]}$, since then both would be $\Rightarrow \lambda = \lambda'$). Since two constructors irrevocably induce different partners for position $i$, their associated structures $S := \phi(T)$ and $S' := \phi(T')$ differ by at least one base
pair, and we have $S \neq S'$, implying the non-ambiguity of $\Delta_{BP}$.

The completeness of $\Delta_{BP}$ requires that any structure $S \in S_n$ can be generated by the evaluation of some term in $\mathcal{T}_q$. This can be established by induction of the length $n \ge \Theta + 2$ of the interval, assuming that, for all $[i', j'] \in Q$ such that $n' := j' - i' + 1 < n$, one has $\Omega_{n'} = S_{n'}$. Now, consider an interval $[i, j]$, $j - i + 1 = n$ and a structure $S \in S_n$, and discuss the partner of $i$ in $S$: if $i$ is unpaired, then $S = \bullet S'$ with $|S'| < n$, so $S'$ is generated by $T' \in \mathcal{T}_{[i+1,j]}$ and $T := \lambda_{\bullet}^{[i,j]}(T') \in \mathcal{T}_{[i,j]}$ such that $\phi(T) = S$; if $i$ is paired to some $k$, then $S = (S') S''$ with $|S'| + |S''| < n$, so $S'$ and $S''$ are generated by terms $T' \in \mathcal{T}_{[i+1,k-1]}$ and $T'' \in \mathcal{T}_{[k+1,j]}$ respectively, so $T := \lambda_k^{[i,j]}(T', T'') \in \mathcal{T}_{[i,j]}$ such that $\phi(T) = S$.

The correctness of $\Delta_{BP}$ requires that, for any secondary structure in the search space, the number of base-pairs is obtained by adding numerical values mapped to constructors. This is possible since, for any term $T$ evaluated as $\phi(T) = S$, the number of base pairs in $S$ coincides with the number of constructors of type $\lambda_k^{[i,j]}$ in $T$, so the scoring function $f: \lambda \to \mathbb{R}$ defined as

$$f(\lambda_{\bullet}^{[i,j]}) = f(\lambda_{\emptyset}^{[i,j]}) = 0 \quad \text{and} \quad f(\lambda_k^{[i,j]}) = \begin{cases} -1 & \text{if } \{w_i, w_k\} \in \mathcal{B} \text{ (valid base pair)} \\ +\infty & \text{otherwise} \end{cases}$$

**Remark 2.1.1:** Note that our assumption of scoring schemes that are additively-defined on constructors/transitions represents a limitations in expressivity in comparison to the more general **evaluation algebra** considered in **algebraic dynamic programming** and its extensions [80, 81, 169, 211]. However, general algebras do not allow a smooth transition from optimization to ensemble analyses, so we (slightly) limit the scope of our framework rather than burden our proofs and theorems. Moreover, the current framework captures, without any complexity overhead, every applications of dynamic programming known to this author in Bioinformatics.

### 2.1.3 Classic optimization

With this final notion of correctness being defined, we can finally turn to a more algorithmic dimension of dynamic programming, initially focusing on optimization problems, an historical focus of dynamic programming since its initial pioneering by Bellman [14].

**Problem 1 (DP-based optimization):**

**Input:** A dynamic programming scheme $\Delta$ and a scoring function $f$, such that $(\Delta, f)$ is correct with respect to an **objective function** $F: \Omega_{\Delta} \to \mathbb{R}$

**Output:** Some element $s^* \in S_{\Delta}$ such that $F(s^*) = \max_{s \in S_{\Delta}} F(s)$

Unsurprisingly, this problem can be solved using dynamic programming, using an algorithm consisting of the following steps:
1. Matrix filling: For all state $q \in Q$, traversed in preorder, compute

$$m_q := \max_{\substack{\lambda \\ q \to q_1, \dots, q_k \in \delta}} f(\lambda) + \sum_{i=1}^{k} m_{q_i} \quad (2.1)$$

2. Backtracking: Return $\mathfrak{B}(q_w)$, recursively defined for any $q \in Q$ as:

$$\mathfrak{B}(q) := \lambda(\mathfrak{B}(q_1), \dots, \mathfrak{B}(q_k)), \text{ if } m_q = f(\lambda) + \sum_{i=1}^{k} m_{q_i}, \quad q \to q_1, \dots, q_k \in \delta$$

Note that the preorder in Step 1 always exists due to the acyclicity of derivations. The time complexity of this step is in $O(|Q| + \alpha^* \times |\delta|)$, $\alpha^* := \max_{d \in \delta} \alpha(d)$ for $\alpha(d)$ the arity of a derivation d, i.e. the number of states in its production. The memoization of computed values $m_q$ requires $\Theta(|Q|)$ memory. These complexities hold for the whole algorithm, since Step 2 involves recursing over at most $|Q|$ states (due to acyclicity), and its complexity is typically orders of magnitude below the requirements of Step 1.

### Example 4: RNA BP folding – Energy minimization

Let us illustrate Equation (2.1) in the context of $\Delta_{BP}$, the DP scheme of RNA BP folding. For base-pair maximization, the **objective function** is $F(S) := |S|$, achieved by a scoring function $f$ such that $f(\lambda_k^{[i,j]}) = \{1 \text{ if } (w_i, w_k) \in \mathcal{B}; -\infty \text{ otherwise}\}$, and $f(\lambda_\bullet^{[i,j]}) = f(\lambda_\emptyset^{[i,j]}) = 0$. We get

$$ \begin{aligned} m_{[i,j]} &= \max \left\{ \begin{array}{ll} f(\lambda_{\emptyset}^{[i,j]}) & \text{if } i + \theta \ge j \quad \triangleright [i, j] \xrightarrow{\lambda_{\emptyset}^{[i,j]}} \emptyset \\ f(\lambda_{\bullet}^{[i,j]}) + m_{[i+1,j]} & \text{if } i + \theta < j \quad \triangleright [i, j] \xrightarrow{\lambda_{\bullet}^{[i,j]}} [i+1, j] \\ \max_k f(\lambda_k^{[i,j]}) + m_{[i+1,j]} + m_{[k+1,j]} & \text{if } i + \theta < j \quad \triangleright [i, j] \xrightarrow{\lambda_k^{[i,j]}} [i+1, k-1], [k+1, j] \end{array} \right. \\ = \max \left\{ \begin{array}{ll} 0 & \text{if } i + \theta \ge j \\ m_{[i+1,j]} & \text{if } i + \theta < j \\ \max_{k=i+\theta+1}^{j} 1 + m_{[i+1,j]} + m_{[k+1,j]} & \text{if } i + \theta < j \land (w_i, w_k) \in \mathcal{B} \end{array} \right. \end{aligned} $$

in which one recognizes the classic DP equation reminded in Section 1.2.3.

## 2.2 Exact computation of Ensemble properties

In many applications of ensemble dynamic programming, one attempts to analyze a specific subset of the search space. Examples abound in RNA bioinformatics where an integer-valued feature function, additively defined with respect to the dynamic programming scheme, partition of the secondary structures with respect to their free-energy, base-pair distance to one or several references structure(s). . .
### 2.2.1 Computing the partition function.

A ubiquitous quantity of interest is the **partition function**, whose definition requires integrating over the whole search space, and is used as a normalization term within ensemble studies. As observed by McCaskill [130], the optimization algorithm of any unambiguous DP scheme can be adapted to compute the partition function, through a simple algebraic substitution. Namely, it suffices to substitute $(min/max, +) \to (+, \times)$, coupled with a suitable exponentiation of energy contributions to compute the partition function from the MFE recursions. This observation has given rise to systematic studies decorrelating the DP scheme from its algebra [80, 97] with a specific focus on semi-ring algebras [135, 138].

**Problem 2 (Partition function):**

**Input:** An unambiguous DP scheme $\Delta$ and a scoring $f$, such that $(\Delta, f)$ is correct with respect to an **energy function** $E: \Omega_{\Delta} \to \mathbb{R}$; $\beta \in \mathbb{R}$ a constant

**Output:** The **partition function** $Z_{\Delta}$ of $\Omega_{\Delta}$, defined as

$$Z_{\Delta} = \sum_{s \in \Omega_{\Delta}} e^{-\beta \cdot E(s)}$$

The above problem can be solved by returning $Z_{\Delta} := Z_{q_w}$, following the computation, for all state $q \in Q$, of

$$Z_q = \sum_{s \in \Omega_q} e^{-\beta \cdot E(s)}.$$

Those quantities can be computed recursively (in preorder), using

$$Z_q := \sum_{\substack{\lambda \\ q \to q_1, \dots, q_k \in \delta}} e^{-\beta \cdot E(\lambda)} \times \prod_{i=1}^{k} Z_{q_i} \quad (2.2)$$

The time complexity of this computation is in $O(|Q| + \alpha^* \times |\delta|)$, $\alpha^*$ being the max arity of a constructor, and requires $\Theta(|Q|)$ memory.

**Example 5: RNA BP folding – Partition function**

A reasonable **energy function** is defined as $E(S) := -|S|$, and implicitly used in Nussinov-Jacobson [143] scheme. It is achieved by an eponymous scoring function $E$ such that $E(\lambda_k^{[i,j]}) = \begin{cases} -1 & \text{if } (w_i, w_k) \in B; +\infty \\ 0 & \text{otherwise} \end{cases}$, and $E(\lambda_{\emptyset}^{[i,j]}) = f(\lambda_{\emptyset}^{[i,j]}) = 0$. We get

$$Z_{[i,j]} := \sum \begin{cases} 1 & \text{if } i + \theta \ge j \\ Z_{[i+1,j]} & \text{if } i + \theta < j \\ \sum_{k=i+\theta+1}^{j} e^{\beta \times Z_{[i+1,j]} \times Z_{[k+1,j]}} & \text{if } i + \theta < j \land (w_i, w_k) \in B. \end{cases}$$
The validity of the final result, i.e. the fact that $Z_{[1,n]} = \sum_{S \in S_n} e^{-\beta \cdot E(S)}$, is a direct consequence of the unambiguity, completeness and correctness properties of $\Delta_{BP}$. The time and space complexities are in $\Theta(n^3)$ and $\Theta(n^2)$ respectively.

### 2.2.2 Probabilities in Boltzmann-Gibbs distributions (inside-outside)

In many relevant contexts, the elements of a search space $\Omega$ can be assumed to follow a **Boltzmann-Gibbs distribution**, where the probability of any $s \in \Omega$ is such that

$$ \mathbb{P}(s) = \frac{e^{-\beta \cdot E(s)}}{\mathcal{Z}} \qquad (2.3) $$

where E represents an energy score, and $\beta$ a constant (analogous to a temperature), and $\mathcal{Z}$ is the partition function. In such a context, the probabilities of individual elements induce average properties (e.g. base-pair probabilities) that are extremely relevant to ensemble analyses.

**Definition 2.2.1 (Unicity of a constructor):** A constructor $\lambda \in \Lambda$ is **unique** within a DP scheme $\Delta$ if it occurs at most once in each term $T \in \mathcal{T}_{\Delta}$.

Under the unicity condition, the probability of a constructor $\lambda$ can be obtained by summing the individual probabilities of search space elements that result from its application (or, equivalently, the terms that contain $\lambda$). Fortunately, this property can be computed efficiently using a suitable dynamic programming scheme, as stated below.

**Problem 3 (Boltzmann probability of constructor(s)):**

**Input:** Unambiguous DP scheme $\Delta$ + scoring $f$, correct w.r.t. energy function E; $\beta \in \mathbb{R}$ a constant; and a **set** $C \subset \Lambda$ **of unique constructors**

**Output:** The Boltzmann probabilities of constructors in C:

$$ \forall \lambda \in C : \mathbb{P}(\lambda \in T) = \frac{\sum_{T \in \mathcal{T}_{\Delta}} e^{-\beta \cdot E(\phi(T))}}{\mathcal{Z}_{\Delta}} \qquad (2.4) $$

This problem, which generalizes the computation of production probabilities in probabilistic context-free grammars, is tackled by a variant of the **inside-outside algorithm** [10, 113]. The algorithm is based on the observation that, for any monitored constructor $\lambda^*$, any term $T^* \in \mathcal{T}_{\lambda^*} := \{T \in \mathcal{T}_{\Delta} | \lambda \in T\}$ can be decomposed into:

1. a derivation $d := q \xrightarrow{\lambda^*} q_1, \dots, q_k$ labeled by an occurrence of $\lambda^*$;

2. an **outside** part, a partial term in $T_q \in \mathcal{T}_{\Delta}$, **truncated** on an occurrence of $q$ (left underived). Let us denote by $\mathcal{D}_q$ the set of outside parts leading to $q$;
3. several **inside** parts, i.e. individual **continuations** of the derivation process, starting from $q_1, \dots, q_k$. Denote as $\mathfrak{J}_{q_1, \dots, q_k}$ the set of inside parts generated from $q_1, \dots, q_k$ respectively;

Under the unicity condition, this decomposition is unambiguous. Moreover, the respective energy contributions of the three parts are independently contributing to the energy, and it follows that $\mathcal{Z}_{\lambda^*}$, the partition function restricted to $\mathcal{T}_{\lambda^*}$, obeys

$$ 
\begin{aligned}
\mathcal{Z}_{\lambda^*} &:= \sum_{T \in \mathcal{T}_{\lambda^*}} e^{-\beta \cdot E(T)} \\
&= \sum_{q \xrightarrow{\lambda^*} q_1, \dots, q_k \in \delta} e^{-\beta \cdot E(\lambda^*)} \times \left( \sum_{T_q \in \mathcal{C}_q} e^{-\beta \cdot E(T_q)} \right) \times \prod_{i=1}^{k} \left( \sum_{T_i \in \mathfrak{J}_{q_i}} e^{-\beta \cdot E(T_i)} \right) \\
&= \sum_{q \xrightarrow{\lambda^*} q_1, \dots, q_k \in \delta} e^{-\beta \cdot E(\lambda^*)} \times \mathcal{Y}_q \times \prod_{i=1}^{k} \mathcal{Z}_{q_i}
\end{aligned} 
\quad (2.5) $$

where $\mathcal{Y}_q := \sum_{T_q \in \mathcal{C}_q} e^{-\beta \cdot E(T_q)}$ is the **outside partition function**. Note that, when $\mathcal{Y}_q$ is known, the above equation allows to simultaneously compute the partition functions $\mathcal{Z}_{\lambda}$ for all monitored unique constructors, through a **single pass over the derivations**.

The computation of $\mathcal{Y}_q$ itself can also be performed by **inverting the dynamic programming scheme**, going from a given state back to the root while allowing the further derivation of siblings found along the way. The outside partition function $\mathcal{Y}_q$ of a node can be computed using dynamic programming using infix order, i.e. starting from the root $q_w$ and processing the ancestors of a node before itself, through

$$ \mathcal{Y}_q := \begin{cases} 1 & \text{if } q = q_w \text{ (root)} \\ \sum_{\substack{q_p \xrightarrow{\lambda} q \in \delta \\ \text{s.t. } q \in q}} e^{-\beta \cdot E(\lambda)} \times \mathcal{Y}_{q_p} \times \prod_{\substack{q' \in q \\ q' \neq q}}^{k} \mathcal{Z}_{q'} & \text{otherwise.} \end{cases} \quad (2.6) $$

Overall, the **inside-outside algorithm** solving Problem 3 can be stated as:

* Using Equation (2.2), compute the inside partition function $\mathcal{Z}_q$ for all state $q \in Q$ in preorder; $\rightarrow \mathcal{O}(|Q| + \alpha^* \times |\delta|)$ time

* Using Equation (2.6), compute the outside partition function $\mathcal{Y}_q$ for all state $q \in Q$ in infix order; $\rightarrow \mathcal{O}(|Q| + \alpha^* \times |\delta|)$ time

* Iterate over derivations $q \xrightarrow{\lambda} q_1 \cdots q_k \in \delta$ to compute $(\mathcal{Z}_{\lambda})_{\lambda \in C}$, initially set to 0. If $\lambda \in C$, update $\mathcal{Z}_{\lambda} \leftarrow \mathcal{Z}_{\lambda} + e^{-\beta \cdot E(\lambda^*)} \mathcal{Y}_q \prod_{i=1}^{k} \mathcal{Z}_{q_i}$; $\rightarrow \mathcal{O}(\alpha^* \times |\delta| + |C|)$ time

* Finally, the algorithm returns $\mathbb{P}(\lambda \in T) := \mathcal{Z}_{\lambda}/\mathcal{Z}_{q_w}, \forall \lambda \in C$.

The algorithm runs in time $\mathcal{O}(|Q| + |C| + \alpha^* \times |\delta|)$, $\alpha^*$ being the maximum arity of a constructor, and requires storage for $\mathcal{O}(|Q| + |C|)$ numbers.
### Example 6: RNA BP folding – BP probabilities

In the context of RNA, the inside/outside algorithm can be used to compute **base-pair probabilities**, as done by McCaskill [130]. Inside contributions/partition functions are computed as detailed in Example 5, using $E(\lambda_k^{[i,j]}) = \{-1 \text{ if } (w_i, w_k) \in \mathcal{B}; +\infty \text{ otherwise}\}$ and $E(\lambda_\emptyset^{[i,j]}) = f(\lambda_\emptyset^{[i,j]}) = 0$. The monitored constructors are all $\lambda_k^{[i,j]}$, each **unique** as a position i cannot be assigned twice.

The outside contributions $Y_q$ are then computed through a specialization of Equation (2.6):

$$ Y_{[i,j]} := \sum_{\substack{i,j \\ \sum_{i'=1}^{n} e^{\beta} \times Y_{[i',j]} \times Z_{[i'+1,i-2]} \\ \sum_{j'=j+1}^{n} e^{\beta} \times Y_{[i,j']} \times Z_{[i+1,j']}}} \begin{cases} 1 & \text{if } [i,j] = [1,n] \\ Y_{[i-1,j]} & \text{if } i-1 \ge 1 \\ \sum_{i'=1}^{n} e^{\beta} \times Y_{[i',j]} \times Z_{[i'+1,i-2]} & \text{if } (w_{i'}, w_{i-1}) \in \mathcal{B} \\ \sum_{j'=j+1}^{n} e^{\beta} \times Y_{[i,j']} \times Z_{[i+1,j']}, & \text{if } (w_{i-1}, w_{j+1}) \in \mathcal{B}. \end{cases} $$

We finally obtain the probabilities of constructors through a specialization of Equation (2.5):

$$ \mathbb{P}(\lambda_k^{[i,j]} \in \mathcal{T}) := \begin{cases} \frac{e^{\beta} \times Y_{[i,j]} \times Z_{[i+1,k-1]} \times Z_{[k+1,j]}}{Z_{[1,n]}} & \text{if } (w_i, w_k) \in \mathcal{B} \\ 0 & \text{otherwise} \end{cases} $$

Summing over all values of j, we get the probability of a base-pairs (i, k), $k - i > 0$

$$ \mathbb{P}((i, k) \in \mathcal{S}) := \sum_{j \ge k}^{n} \mathbb{P}(\lambda_k^{[i,j]} \in \mathcal{T}) $$

The probabilities of all base pairs can then be computed the inside/outside, in $\Theta(n^3)$ time, and $\Theta(n^2)$ space by computing the probabilities of constructors on the fly within the above sum.

**Remark 2.2.1 (Beyond unique features):** Remark that, in the case where $\lambda$ is not unique, i.e. it occurs more than once in a term, the output of the above algorithm is no longer the probability of occurrence, but the **expected number of occurrences** of $\lambda$ in a random, Boltzmann-distributed, term. This quantity may be of interest, for instance when trying to assess expected properties of the Bolzmann ensembles, since it allows the simultaneous computation of many expected features in a single pass.

### 2.2.3 Ensemble centroid and Maximum Expected Accuracy (MEA)

The probabilities computed in Problem 3 allow to assess a notion of **support** for the individual **features** (e.g. base pairs, helices...) of a solution within the Boltzmann-Gibbs distribution. Thus, they can be used to assess how **representative** a given solution is of the Boltzmann-Gibbs ensemble. In particular, when $s^*$ is the minimum free-energy solution, we know that $s^*$ achieves maximal probability in the Boltzmann-Gibbs distribution.

However, in absolute terms, the probability of $s^*$ may be (and usually is) abysmally small, and does not allow in itself to distinguish between two very different situations:

1. The solution $s^*$ is surrounded by a family of **similar suboptimal solutions** {$s'_1, s'_2, \dots$}, having very similar features (e.g. $|s^*, s'_i| \le \eta$ for some notion of distance), **overtak-**
ing the probability distribution ($\mathbb{P}(s^*) + \sum_i \mathbb{P}(s'_i) = 1 - \varepsilon$);

2. The solution $s^*$, even supplemented by similar suboptimals, is **highly dominated** by dissimilar solutions in the Boltzmann ensemble ($\mathbb{P}(s^*) + \sum_i \mathbb{P}(s'_i) = \varepsilon$).

A possible way to distinguish between those two worlds, consists in computing an **expected distance** of $s^*$ to a random solution in the ensemble.

**Definition 2.2.2 (Weighted distance between solutions):** Given two solutions $s, s' \in \Omega$ resulting from the applications of sets of unique constructors $C = \{\lambda_1, \dots, \lambda_k\}$ and $C' = \{\lambda'_1, \dots, \lambda'_{k'}\}$ respectively. Then the **distance** $|s, s'|$ **between** $s$ and $s'$ is defined as

$$|s, s'|_{\mu} = \sum_{\lambda \in \Lambda} \pi_{\lambda} \times (\mathbf{1}_{\lambda \in C} - \mathbf{1}_{\lambda \in C'})^2$$

where $\pi: \Lambda \to \mathbb{R}^+$ is a collection of weights.

Since constructors represent atomic operations that build a given element of the search space (e.g. adding a base pair, declaring a nucleotide unpaired...), this notion of distance represents a natural way to represent popular distance metrics.

Equipped with a notion of distance, we can now define the **centroid** of the Boltzmann-Gibbs ensemble [55, 87] as its most central element, i.e. the solution having **minimum** **expected distance** to a, Boltzmann-distributed, elements of the search space.

**Problem 4 (Centroid solution):**

**Input:** Unambiguous DP scheme $\Delta$ + scoring $f$, correct w.r.t. function $E$; $\beta \in \mathbb{R}$ a constant; and a weighted distance $|\star, \star|_{\mu}$

**Output:** Solution $s^* \in \Omega_{\Delta}$ minimizing the **expected distance to the ensemble**:

$$s^* = \underset{s \in \Omega_{\Delta}}{\operatorname{argmin}} \sum_{s \in \Omega_{\Delta}} \mathbb{P}(s') \times |s, s'|_{\mu} \quad (2.7)$$

Fortunately, the expected distance to the ensemble of any given candidate solution $s \in \Omega$ can be reexpressed as a simple sum over the Boltzmann probabilities of constructors.
Indeed, one has

$$
\begin{align*}
\sum_{s' \in \Omega_{\Delta}} \mathbb{P}(s') \times |s, s'|_{\mu} &= \sum_{s' \in \Omega_{\Delta}} \mathbb{P}(s') \times \sum_{\lambda \in \Lambda} \pi_{\lambda} \times (\mathbf{1}_{\lambda \in s} - \mathbf{1}_{\lambda \in s'})^2 \\
&= \sum_{\substack{\lambda \in \Lambda \\ \lambda \in s}} \pi_{\lambda} \times \mathbb{P}(\lambda \notin T) + \sum_{\substack{\lambda \in \Lambda \\ \lambda \notin s}} \pi_{\lambda} \times \mathbb{P}(\lambda \in T) \\
&= \sum_{\substack{\lambda \in \Lambda \\ \lambda \in s}} \pi_{\lambda} \times (\mathbb{P}(\lambda \notin T) - \mathbb{P}(\lambda \in T)) + \sum_{\lambda \in \Lambda} \pi_{\lambda} \times \mathbb{P}(\lambda \in T) \quad (2.8)
\end{align*}
$$

Since the rightmost sum no longer depends on s, finding the solution that minimizes
the expected distance is equivalent to finding a solution that optimizes the leftmost sum.
Following this observation, one can solve Problem 4 by executing the following algorithm:

1. Compute the Boltzmann probabilities of constructors used by the weighted distance, i.e. solve Problem 3 with $C := \{\lambda \in \Lambda \mid \pi_{\lambda} \neq 0\}$;

2. Find $s \in \Omega_{\Delta}$ that minimizes the leftmost term of Equation 2.8, i.e. solve Problem 1, maximizing the objective function $F: \Lambda \to \mathbb{R}$ such that

$$
F(\lambda) := \pi_{\lambda} \times (\mathbb{P}(\lambda \in T) - \mathbb{P}(\lambda \notin T)) = \pi_{\lambda} \times (2 \mathbb{P}(\lambda \in T) - 1) \quad (2.9)
$$

Example 7: RNA BP folding – Centroid computation

We consider the classic **base-pair distance** as the distance to be minimized, and accordingly replace all constructors $\lambda_k^{[i,j]}$ with new **simplified constructors** $\lambda_{(i,k)}$ and $\lambda_{(i)}$ which respectively represent occurrences of a base pair (i, k) and an unpaired position i, irrespectively of their context [i, j] of creation (since the context of a base pair should not contribute to the distance). We set the weight of all constructors to 0 except for $\pi(\lambda_{(i,k)}) = 1$ in the distance definition, and compute the base-pair probabilities $p_{i,j} := \mathbb{P}((i, j) \in S)$ as shown in Section 2.2.2.

Then, we solve Problem 1 in this new setting, i.e. compute the recurrence

$$
c_{[i,j]} = \max \begin{cases} 0 & \text{if } i + \theta \ge j \\ c_{[i+1,j]} & \text{if } i + \theta < j \\ \max_{k=i+\theta+1}^{j} (2 p_{i,k} - 1) + c_{[i+1,j]} + c_{[k+1,j]} & \text{if } i + \theta < j \land (w_i, w_k) \in \mathcal{B} \end{cases}
$$

to get the least distance of a structure to the ensemble (up to a constant, i.e. the rightmost term in (2.9))
A classic backtrack allows to recover the **centroid secondary structure**.

A **Maximum Expected Accuracy (MEA)** solution [121] can be obtained in a very similar fashion by simplifying the objective function of Equation (2.9) to $F(\lambda) := \pi_{\lambda} \times \mathbb{P}(\lambda \in T)$, with $\pi_{\lambda_{BP}} := 2\gamma$ for base-pairing constructors, and $\pi_{\lambda_{Unp.}} := 1$ for unpaired constructors.

2.2.4 General moments of additive scores (150)

Given an additive scoring function F, it is a natural question to ask for the induced
distribution of F under a Boltzmann Gibbs distribution. Since most such distributions
are typically Gaussian, a first task is to compute the expected value μ<sub>Δ</sub>(F) and variance
$\sigma_{\Delta}(F)$ of $F$, respectively defined as $\mu_{\Delta}(F) := \mathbb{E}(F(T))$ and

$$ \sigma_{\Delta}(F) := \sqrt{\sum_{T \in \mathcal{T}_{\Delta}} \mathbb{P}(T) \cdot (F(T) - \mu_{\Delta}(F))^2} = \sqrt{\mathbb{E}(F(T)^2) - \mathbb{E}(F(T))^2}. $$

In order to capture characteristics of more general distributions, one may consider the **moments** of the distribution, defined as $\mathbb{E}(F(T))$, $\mathbb{E}(F(T^2))$, ..., previously considered by Miklos, Meyer and Borbala [136] specifically for the free-energy. The **cross-moments** $\mathbb{E}(F_1(T)^{n_1} \cdot F_2(T)^{n_2} \cdots)$ of multiple functions are also of potential interest, as their computation allows to derive the **Pearson correlation** $\rho_{\Delta}(F_1, F_2)$ of two functions $F_1$ and $F_2$ through

$$ \begin{aligned} \rho_{\Delta}(F_1, F_2) &= \frac{\mathbb{E}((F_1(T) - \mu_{\Delta}(F_1)) \times (F_2(T) - \mu_{\Delta}(F_2)))}{\sigma_{\Delta}(F_1) \times \sigma_{\Delta}(F_2)} \\ &= \frac{\mathbb{E}(F_1(T) \times F_2(T)) - \mathbb{E}(F_1(T)) \times \mathbb{E}(F_2(T))}{\sqrt{\mathbb{E}(F_1(T)^2) - \mathbb{E}(F_1(T))^2} \times \sqrt{\mathbb{E}(F_2(T)^2) - \mathbb{E}(F_2(T))^2}} \end{aligned} $$

so the correlation can be computed from the evaluation of $\mathbb{E}(F_1(T)^{n_1} \times F_2(T)^{n_2})$ for all values of $(n_1, n_2) \in \{(1, 0), (0, 1), (2, 0), (0, 2), (1, 1)\}$.

In a collaboration with Cédric Saule [150], we have considered the computation of general cross-moments within dynamic programming schemes.

**Problem 5 (Cross moments of a DP scheme):**

**Input:** Unambiguous DP scheme $\Delta + \text{scoring } f$, correct w.r.t. function $E$; $\beta \in \mathbb{R}$; and **scoring functions** $(F_1, \dots, F_p)$ with associated **degrees** $(\gamma_1^*, \dots, \gamma_p^*)$

**Output:** The (cross) moment $\gamma_1^*, \dots, \gamma_p^*$ for $F_1, \dots, F_p$:

$$ m^{\{\gamma_1^* \cdots \gamma_p^*\}} := \mathbb{E} \left( F_1(T)^{\gamma_1^*} \times \cdots \times F_p(T)^{\gamma_p^*} \right) = \sum_{T \in \mathcal{T}_{\Delta}} \mathbb{P}(T) \prod_{i=1}^{p} F_i(T)^{\gamma_i^*} \quad (2.10) $$

A reexpression of our result [150] states that $m^{\gamma_1, \dots, \gamma_p}$ can be computed as

$$ m_{q}^{\{\gamma_{1}^{*} \cdots \gamma_{p}^{*}\}} := \sum_{\lambda \to q_{1} \cdots q_{k} \in \delta} e^{-\beta \cdot E(\lambda)} \sum_{\substack{\lambda_{1} + \tau_{1,1} \cdots \tau_{1,k} = \gamma_{1} \\ \vdots \\ \lambda_{p} + \tau_{p,1} \cdots \tau_{p,k} = \gamma_{p}}} \prod_{i=1}^{p} F_{i}(\lambda)^{\lambda_{i}} \left( \prod_{\tau_{i,1} \cdots \tau_{i,k}} \gamma_{i} \right) \prod_{j=1}^{k} m_{q_{j}}^{\{\tau_{1,j} \cdots \tau_{p,j}\}}. \quad (2.11) $$

This expression, whose underlying intuition is probably hard to decipher at first sight, is in fact strongly inspired by the **partial pointing operator** in enumerative combinatorics [50, 69]. Its main underlying idea is to modify the DP scheme $\Delta$, introducing of a **controlled ambiguity**, into a DP scheme $\Delta^{\{\gamma_1, \dots, \gamma_p\}}$ designed to generate pointed/weighted versions of the original terms.
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>Figure 2.1: Distribution of base pairs induced by a uniform (left) and Boltzmann distribution. Empirical distributions (histograms) measured on 20 000 stochastic backtracks, approximated by a Normal distribution with exact means and variances, computed from the first two moments. Boltzmann Sampling performed at 37°C with E(S) = -|S| (right) for a toy RNA w := (ACGU)²⁵.</figcaption></figure>

By generating each term of Δ with suitable multiplicity, and weighting them with a carefully chosen monomial in $F_1, F_2, \dots, \Delta^{\{\gamma_1, \dots, \gamma_p\}}$ we ensure that the **multinomial formula** applies and, for all $T \in \Delta$, the overall weight contributed by the equivalence class $\{T' \in T_{\Delta^{\{\gamma_1, \dots, \gamma_p\}}} | S' \Leftrightarrow S\}$ simplifies into $F_1^{\gamma_1}(T) \times \dots \times F_p^{\gamma_p}(T)$.

An algorithm for Problem 5 then consists in using Equation (2.11) to compute $m_q^{\gamma_1 \cdots \gamma_p}$ for all $(\gamma_1 \cdots \gamma_p) \le (\gamma_1^* \cdots \gamma_p^*)$ in increasing lexicographic order and all $q \in Q$ in preorder for each degree vector, and return $m_{q_w}^{\{\gamma_1^* \cdots \gamma_p^*\}}$. Its time complexity is loosely bounded by $O(|\delta| \times \prod_{i=1}^p \gamma_i^{1+\alpha^*} \times P(k, p))$ where P is a polynomial of bounded degree, and the space complexity is in $\Theta(|Q| \times \prod_{i=1}^p \gamma_i^*)$. In practice, degree vectors are typically small values, so the complexity of the algorithm is equivalent, up to a constant factor, to the computation of the partition function.

### Example 8: RNA BP folding – Distribution of base pairs

The homopolymer model predicts that the distribution of base-pairs in random, Boltzmann-distributed secondary structure asymptotically follows a normal distribution, and is thus characterized by its first two moments $m^{\{1\}}/\mathcal{Z}$ and $m^{\{2\}}/\mathcal{Z}$. Those can be computed by introducing a scoring function F, defined as $F(\lambda^{[i,j]}) = F(\lambda_{\emptyset}^{[i,j]}) = 0$ and $F(\lambda_k^{[i,j]}) = 1$, leading to

$$ m_{[i,j]}^{\{1\}} = \begin{cases} F(\lambda_{\emptyset}^{[i,j]}) & \text{if } i + \theta \ge 1 \\ F(\lambda_{\bullet}^{[i,j]}) \times m_{[i,j]}^{\{0\}} + m_{[i,j]}^{\{1\}} & \text{if } i + \theta < 1 \\ \sum_{k=i+\theta+1}^{j} F(\lambda_k^{[i,j]}) \times e^{\beta} \times m_{[i+1,j]}^{\{0\}} \times m_{[k+1,j]}^{\{0\}} & \text{if } i + \theta < j \land (w_i, w_k) \in \mathcal{B} \\ \sum_{k=i+\theta+1}^{j} e^{\beta} \times m_{[i+1,j]}^{\{1\}} \times m_{[k+1,j]}^{\{0\}} & \text{if } i + \theta < j \land (w_i, w_k) \in \mathcal{B} \\ \sum_{k=i+\theta+1}^{j} e^{\beta} \times m_{[i+1,j]}^{\{0\}} \times m_{[k+1,j]}^{\{1\}} & \text{if } i + \theta < j \land (w_i, w_k) \in \mathcal{B} \end{cases} $$

with $m_{[i,j]}^{\{0\}} := \mathcal{Z}_{[i,j]}$, leading to the exact value for the expected number of base pairs $\mu_{BP} :=$
$m_{[1,n]}^{(1)}/\mathcal{Z}_{[1,n]}$. A similar equation can be derived for the second moment

$$m_{[i,j]}^2 = \begin{cases} 0 & \text{if } i + \theta \ge 1 \\ m_{[i,j]}^2 & \text{if } i + \theta < 1 \\ \sum_{k=i+\theta+1}^j e^\beta \times F(\lambda_k^{[i,j]})^2 \times m_{[i+1,j]}^0 \times m_{[k+1,j]}^0 & \text{if } i + \theta < j \land (w_i, w_k) \in \mathcal{B} \\ \sum_{k=i+\theta+1}^j e^\beta \times 2 \times F(\lambda_k^{[i,j]})^1 \times m_{[i+1,j]}^1 \times m_{[k+1,j]}^0 & \text{if } i + \theta < j \land (w_i, w_k) \in \mathcal{B} \\ \sum_{k=i+\theta+1}^j e^\beta \times 2 \times F(\lambda_k^{[i,j]})^1 \times m_{[i+1,j]}^0 \times m_{[k+1,j]}^1 & \text{if } i + \theta < j \land (w_i, w_k) \in \mathcal{B} \\ \sum_{k=i+\theta+1}^j e^\beta \times 2 \times m_{[i+1,j]}^1 \times m_{[k+1,j]}^0 & \text{if } i + \theta < j \land (w_i, w_k) \in \mathcal{B} \\ \sum_{k=i+\theta+1}^j e^\beta \times m_{[i+1,j]}^2 \times m_{[k+1,j]}^0 & \text{if } i + \theta < j \land (w_i, w_k) \in \mathcal{B} \\ \sum_{k=i+\theta+1}^j e^\beta \times m_{[i+1,j]}^0 \times m_{[k+1,j]}^2 & \text{if } i + \theta < j \land (w_i, w_k) \in \mathcal{B} \end{cases}$$

from which we get the standard deviation $\sigma_{BP} := \sqrt{m_{[1,n]}^2/\mathcal{Z}_{[1,n]} - \mu_{BP}^2}$

### 2.2.5 Classified DP with the Discrete Fourier Transform (DFT) (173, 174)

In some situations, precise aspects of the distribution are of interest, and cannot be easily captured by summary statistics. For instance, Freyhult et al [75] partition the Boltzmann Ensemble according to the distance to a reference secondary structure. Considering the distribution of distance, they observe bimodal (an even trimodal) distributions that they interpret as potential evidences for the presence of multistable RNAs. Other examples in RNA bioinformatics include the classification of sequences/structures with respect to their Hamming distance to a wild-type sequence, to assess the mutational robustness of concrete RNAs [191], the computation of the density of states, i.e. the distribution of free-energies within an RNA [43], or the projection of energy landscapes with respect to the unfolded and native states [119].

Those instances can be expressed as a specialized version of **classified dynamic programming**, where additional parameters are used to partition the search space.

**Problem 6 (Classified partition function):**

**Input:** Unambiguous DP scheme $\Delta$ + scoring $f$, correct w.r.t. function $E$; $\beta \in \mathbb{R}$; and an additive **scoring functions** $F: \Lambda \to \mathbb{N}^+$ with $V := \max_{T \in \mathcal{T}} F(T)$

**Output:** The **classified partition function** $\mathcal{Z}^{[v]}$, $0 \le v \le V$, defined as

$$\mathcal{Z}^{[v]} := \mathbb{P}(F(T) = v) \times \mathcal{Z}_{\Delta} = \sum_{\substack{T \in \mathcal{T} \\ \text{s.t. } F(T) = v}} e^{-\beta \cdot E(T)}$$

A classic approach, used in many works including the above references, would consist in adapting the dynamic programming scheme through **explicit convolution products**,
leading to the following recurrence:

$$ \mathcal{Z}_q^{[v]} := \sum_{\substack{\lambda \\ q \to q_1, \dots, q_k \in \delta}} \sum_{v_1 + \dots + v_k + F(\lambda) = v} e^{-\beta \cdot E(\lambda)} \prod_{i=1}^{k} \mathcal{Z}_{q_i}^{[v_i]} \quad (2.12) $$

While conceptually simple, this approach induces extreme computational demands, with a time complexity in $\Theta(|\delta| \cdot V^{\alpha^*})$, $\alpha^*$ being the max arity of a derivation, and using $\Theta(|Q| \cdot V)$ memory. For instance, in our running DP scheme example for RNA secondary structures (see Section 1.2.3), classifying according to the number of base pairs and applying this strategy induces an algorithm in $\Theta(n^5)$ time and $\Theta(n^3)$ memory complexities.

In collaboration with Evan Senter, Ivan Dotu, Peter Clote and Saad Sheikh, we have proposed an alternative strategy based on the Discrete Fourier Transform (DFT) to avoid costly convolution products [173, 174], following observations made in an earlier collaborative work with Jérôme Waldispühl [192, 193].

Its core idea is to consider a generalized version of the partition function, seen as a polynomial in a formal variable $x$, which tracks the increments of the scoring function, such that

$$ \mathcal{Z}(x) = \sum_{T \in T_{\Delta}} e^{-\beta \cdot E(T)} x^{F(T)} = \sum_{v=0}^{V} \mathcal{Z}^{[v]} x^{v}. $$

Such a polynomial can be evaluated at any point without having to previously determine its coefficients, using the following DP equation

$$ \mathcal{Z}_q(x) := \sum_{\substack{\lambda \\ q \to q_1, \dots, q_k \in \delta}} e^{-\beta \cdot E(\lambda)} x^{F(\lambda)} \prod_{i=1}^{k} \mathcal{Z}_{q_i}(x). \quad (2.13) $$

From a well-chosen set of evaluations of $\mathcal{Z}(x) := \mathcal{Z}_{q_w}(x)$, it is then possible to use **polynomial interpolation** to recover the coefficients of the polynomial, i.e. the classified partition functions $\mathcal{Z}^{[v]}$.

For the sake of illustration, one may perform a preliminary evaluation of $\mathcal{Z}$ at $(1+V)$ distinct points $x_0, x_1, \dots, x_V$, and then use **Gaussian elimination** to solve the linear system

$$ \begin{bmatrix} x_0^V & \cdots & x_0^1 & 1 \\ x_1^V & \cdots & x_1^1 & 1 \\ \vdots & & \vdots & \vdots \\ x_V^V & \cdots & x_V^1 & 1 \end{bmatrix} \times \begin{bmatrix} \mathcal{Z}^{[0]} \\ \mathcal{Z}^{[1]} \\ \vdots \\ \mathcal{Z}^{[V]} \end{bmatrix} = \begin{bmatrix} \mathcal{Z}(x_0) \\ \mathcal{Z}(x_1) \\ \vdots \\ \mathcal{Z}(x_V) \end{bmatrix}. \quad (2.14) $$

in its unknowns $\mathcal{Z}^{[0]}, \mathcal{Z}^{[1]}, \dots$ thus solving Problem 6 in $\Theta(V|\delta|\alpha^* + V^3)$ time. However, using Gaussian elimination would be highly impractical due to its many high **numerically instability**. Moreover, the $\Theta(V^3)$ time complexity of elimination alone would become
<caption>**Algorithm 1:** Fast Fourier Transform</caption>

**Result:** The Fourier transform of x, a vector [$y_j = \sum_{k=0}^{m-1} x_k \omega_m^{jk} | j \in [0, m-1]$]

**Function** FFT(Coefficients x, m = |x| = 2^d):

if m = 1 then return y;

$\omega \leftarrow 1;$

($\mathbf{x}^{\text{even}}, \mathbf{x}^{\text{odd}} \leftarrow ([x_0, \dots, x_{m-2}], [x_1, \dots, x_{m-1}]);$

($\mathbf{y}^{\text{even}}, \mathbf{y}^{\text{odd}} \leftarrow (\text{FFT}(\mathbf{x}_{\text{even}}), \text{FFT}(\mathbf{x}_{\text{odd}}));$

for $k \leftarrow 1$ to $m/2 - 1$ do

$\mathbf{y}_k \leftarrow \mathbf{y}_k^{\text{even}} + \omega \times \mathbf{y}_k^{\text{odd}};$

$\mathbf{y}_{k+m/2} \leftarrow \mathbf{y}_k^{\text{even}} - \omega \times \mathbf{y}_k^{\text{odd}};$

$\omega \leftarrow \omega \times \omega_n;$

return [$y_0, y_1, \dots, y_{m-1}$]

unreasonable for larger ranges of values for F, needed to emulate several scoring functions
as pointed out in Remark 2.2.2 below.

Enters the **(Inverse) Discrete Fourier Transform (IDFT)** [41, 42] which allows, in $\mathcal{O}(m \log m)$
time, to recover the coefficients $\mathbf{x} := (x_0, \dots, x_{m-1})$ of a polynomial $P(x)$ of degree $m-1$,
$m = 2^d$, from its evaluation at the $m$-th **roots of the unity**. Indeed, let $\omega_m := e^{2\pi i/m}$,
then the relationship between $\mathbf{x}$ and $\mathbf{y} := \{P(\omega_m^k)\}_{k=0}^{m-1}$ remarkably simplifies

$$ \forall 0 \le k < m : x_k = \frac{1}{m} \sum_{j=0}^{m-1} y_j \times \omega_m^{-j} \quad (2.15) $$

Moreover, evaluating all coefficients in $\mathbf{x}$ can be achieved much faster than the $\Theta(n^2)$
algorithm suggested by the above equation, using the **Fast Fourier Transform** to speed
up the evaluation to $\Theta(n \log n)$.

To use the inverse DFT and interpolate $\mathcal{Z}(\mathbf{x})$, thus efficiently solving Problem 6, we first
trivially extend Equation (2.13) to support complex arguments. Then we proceed to the
main algorithm:

1. Round $V + 1$ upwards to $m = 2^d$, $d \in \mathbb{N}$, the closest power of 2;

2. Evaluate $\mathcal{Z}(\omega_m^0), \mathcal{Z}(\omega_m^1), \dots, \mathcal{Z}(\omega_m^{m-1}) \to y_0, y_1, \dots, y_{m-1}$ with $\omega_m = e^{2\pi i/m}$;

3. Apply the FFT Algorithm 1 to [$y_0, y_1, \dots, y_{m-1}$] to obtain

$$ \left[ x'_{k} := \sum_{j=0}^{m-1} y_{j} \omega_{m}^{j k} \right]_{k=0}^{m-1} . $$

Notice from Equation (2.15) that $x'_{m-k} = m x_k, \forall k \in [1, m-1]$ and $x'_0 = m x_0$;

4. Return $\left[ \mathcal{Z}^{[v]} := \frac{x'_{\phi(v)}}{n} \mid v \in [0, V], \phi(0) = 0 \text{ and } \phi(v > 0) = V + 1 - v \right]$.
This algorithm has time complexity in $Θ(V|δ|α^* + V \log V)$ and space complexity in $Θ(|Q| + V)$. Moreover, being dominated by the independent evaluations of the $Z(ω_m^i)$ terms, it can be largely accelerated by a trivial parallel execution, e.g. leading to a runtime in $Θ(|δ|α^* + V \log V)$ on V processors.

**Remark 2.2.2 (Multidimensional classified DP):** The DFT-based algorithm can easily be extended to support multiple scoring functions $F_1, \dots, F_p$, associated to maximum values $V_1, \dots, V_p$, to ultimately compute $Z^{[v_1, \dots, v_p]}$ the partition function classified by values $(v_1, \dots, v_p)$ of the scoring functions.

This use-case can indeed be emulated using only a **single function** $F$ defined on constructors as $F(λ) = \sum_{i=1}^{p} F_i(λ) \prod_{j=1}^{i-1} (1 + V_j)$, and extended additively into

$$F(T) = \sum_{i=1}^{p} F_i(T) \prod_{j=1}^{i-1} (1 + V_j) = F_1(T) + (1 + V_1) \times F_2(T) + (1 + V_2) \times F_3(T) \dots$$

In this setting, any value $v$ for $F$ **encodes** a vector $(v_1, \dots, v_p)$, such that the value of $Z^{[v_1, \dots, v_p]}$ can be read in $Z^{[v]}$ after the FFT determination of $1 + V = \prod_i (1 + V_i)$ coefficients.

## 2.3 Probabilistic estimates

**Sampling methods** provide an alternative to exact computation, allowing the estimation of statistical ensemble properties for arbitrary features. Namely, general properties of Boltzmann ensembles can also be **estimated** from a **statistically representative sample** of candidate solutions, following the approach introduced by Ding and Lawrence [51] in the context of RNA folding prediction.

Considering a general **feature function** $F: Ω \to R$, one first **generates** a sequence $s = (s_1, s_2, \dots, s_M) \in Ω^*$ of M random elements from the search space, each drawn from the Boltzmann distribution, and then returns the **empirical mean**

$$\hat{F}(s) = \frac{\sum_{s \in s} F(s)}{M} \quad \text{such that} \quad \lim_{|s| \to \infty} \hat{F}(s) = \mathbb{E}(F(S)). \qquad (2.16)$$

For instance, the Boltzmann probability of a base-pair $(i, j)$ can be estimated by considering a **Boolean feature function** $F_{(i,j)}(s) = 1$ if $(i, j) \in s; 0$ otherwise. More generally, this strategy allows to capture feature functions of arbitrary complexity, including those that are not additive with respect to any dynamic programming scheme.

A first algorithmic difficulty lies in the generation itself and is usually tackled through the simulation of a well-calibrated Markov process in classic Bayesian inference [5, 76]. However, complex distributions may induce high mixing times, before the process converges to its steady-state. Moreover, to ensure that the targeted distribution is indeed sampled (i.e. that the simulation has been executed sufficient long to guarantee converge)
requires highly-technical analyses [152] that do not appear to lend themselves to generic algorithmic design.

This section focuses on an alternative, where an unambiguous dynamic programming scheme for the partition function (but not necessarily with respect to the feature function) is used for the generation. Probabilities for each of the transitions are precomputed, so that generated objects are independent and follow the targeted distribution. This can be seen as an instance of the **recursive method** in random generation [71, 197].

### 2.3.1 Foreword: On the number of samples (165)

The choice of the **number of samples** is critical when using sampling to estimate statistical properties. This number should be large enough to yield accurate estimates, but low enough to preserve efficiency. Historically, and in many subsequent works, a sample size of 1 000 structures was proposed [51], somewhat irrespectively of the precise context. However, such a *one size fits all* approach may not yield accurate results, motivating the following discussion and recommendations.

The **empirical mean** estimator represents a sum of independent variables, meaning that classic concentration inequalities apply with minimal modifications. In particular, the **Hoeffding inequality** implies that, for any feature F:

$$P(|\widehat{F}(S) - \mathbb{E}(F(S))| \ge \varepsilon) \le 2 \exp\left(\frac{-2m\varepsilon^2}{c}\right), \quad (2.17)$$

where $\varepsilon$ is a tolerated **absolute error** level, $S$ is a random sample of size $m$ and $c := (\max_S(F(S)) - \min_S(F(S)))^2$ is a trivial upper bound of the variance of the feature. Note that when a feature function takes binary values 0/1, e.g. when estimating a probability, then one has $c = 1$. Equation (2.17) can be used to build a **confidence interval** at level $(1 - \alpha)$, for any value $\alpha \in [0, 1]$, and we get:

$$\left[ \widehat{F}(S) - \sqrt{\frac{c}{2m} \log \left( \frac{2}{\alpha} \right)}, \widehat{F}(S) + \sqrt{\frac{c}{2m} \log \left( \frac{2}{\alpha} \right)} \right].$$

This means that, over multiple estimations from sampled elements, at least a fraction $(1 - \alpha)$ of the runs will produce errors smaller than $\sqrt{\frac{c}{2m} \log \left( \frac{2}{\alpha} \right)}$. This function can be inverted (numerically) to estimate the number $m$ of samples that achieve an **absolute error** bounded by $\varepsilon$ at least $(1 - \alpha)$ of the times.

We report in Table 2.1 typical sample sizes required to achieve a given precision with reasonable probability when estimating probabilities (i.e. expectations of 0/1-valued features). For instance, to reach a **90% chance** of estimating a base pair probability **within 0.5% of its true value**, a total of **59 915 structures** should be generated. In
<table><caption>Table 2.1: **Recommended number of samples to estimate probabilities (boolean features).** For instance, to ensure that the estimate falls within 1% of the true value for 95% of the runs, a large number of m = 18 444 structures should be generated.</caption><thead><tr><th rowspan="2">Tolerated Error</th><th colspan="3">Frequency within tolerance</th></tr><tr><th>90%</th><th>95%</th><th>99%</th></tr></thead><tbody><tr><td>ε = 20%</td><td>37</td><td>46</td><td>66</td></tr><tr><td>ε = 10%</td><td>150</td><td>184</td><td>265</td></tr><tr><td>ε = 5%</td><td>599</td><td>738</td><td>1 060</td></tr><tr><td>ε = 2.5%</td><td>2 397</td><td>2 951</td><td>4 239</td></tr><tr><td>ε = 1%</td><td>14 979</td><td>18 444</td><td>26 492</td></tr><tr><td>ε = 5‰</td><td>59 915</td><td>73 778</td><td>105 966</td></tr><tr><td>ε = 1‰</td><td>1 497 866</td><td>1 844 440</td><td>2 649 159</td></tr></tbody></table>

In particular, **1 000 structures** norm, usually considered in the literature, will guarantee a value **within 3%** of the true probability **only 2/3 of the times**, although this sample size will **almost always** (99%) return estimates **within 5%** of the correct value.

### 2.3.2 Statistical sampling

Let us now turn to the problem of sampling from the Bolzmann-Gibbs distribution, which we formalize as follows.

**Problem 7 (Statistical sampling):**

**Input:** Unambiguous DP scheme $\Delta$ + scoring $f$, correct w.r.t. function $E$; $\beta \in \mathbb{R}$ a constant; $M$ a **number of samples**

**Output:** A sequence $s \in \Omega_{\Delta}^{*}$, $|s| = M$ of elements of the search space, independently generated in the Boltzmann distribution

$$P(s) = \frac{e^{-\beta \cdot E(s)}}{Z_{\Delta}}, \forall s \in s.$$

In the context of an unambiguous/complete DP scheme, one can easily adapt the partition function recurrences into a **stochastic backtrack algorithm** to solve the statistical sampling problem, as done by Ding and Lawrence [51] in the context of RNA. First, the algorithm performs a preliminary computation of the partition function $Z_q$ for all states $q \in Q$, using the DP algorithm introduced for Problem 2. The algorithm then generates $M$ independent random terms/solutions using a **stochastic backtrack** procedure $\tilde{B}$ described in Algorithm 2.

The correctness of the algorithm can be established by a simple induction. Indeed, consider a state $q^* \in Q$ and assume that, invoked on any state accessible from $q$, the
<figure><img src="image_19.png" /><figcaption>**Algorithm 2:** Stochastic backtrack</figcaption></figure>

algorithm generates a term $T \in \mathcal{T}_q$ with probability

$$ \mathbb{P}(T | q) = \frac{e^{-\beta \cdot E(T)}}{Z_q} $$

Executing Algorithm 2 on $q^*$, a derivation $d := q^* \xrightarrow{\lambda} q_1, \dots, q_k$ is chosen with probability $e^{-\beta \cdot E(\lambda)} \prod_{i=1}^k Z_{q_i}/Z_{q^*}$, followed by the independent recursive generation of terms $T_1, \dots, T_k$ from each of the states $q_1, \dots, q_k$. Due to the unambiguity of $\Delta$, the resulting term $T^* := \lambda(T_1, \dots, T_k)$ cannot be generated from any other derivation of $q^*$, and its **emission probability** is thus

$$ \begin{align*} \mathbb{P}(T^* | q^*) &= \frac{e^{-\beta \cdot E(\lambda)} \prod_{i=1}^k Z_{q_i}}{Z_{q^*}} \times \frac{e^{-\beta \cdot E(T_1)}}{Z_{q_1}} \times \dots \times \frac{e^{-\beta \cdot E(T_k)}}{Z_{q_k}} \\ &= \frac{e^{-\beta \cdot (E(\lambda) + E(T_1) + \dots + E(T_k))}}{Z_{q^*}} \\ &= \frac{e^{-\beta \cdot E(\lambda(T_1, \dots, T_k))}}{Z_{q^*}} \end{align*} $$

in which one recognizes the targeted distribution. The induction is completed by an inspection of leaf states, whose derivations do not require further recursive calls, to verify that they indeed locally induce a Boltzmann distribution.

The time complexity of the algorithm can be trivially bounded by $O(M \cdot |\delta| \cdot \alpha^*)$, with $\alpha^*$ being the max arity of a derivation, but is usually substantially lower since this upper bound would unrealistically require all states to be traversed by all backtracks. A more refined analysis could introduce upper bounds respectively on the numbers of states traversed during a generation ($c^+$), and derivations available from a state ($d^+$). Since one typically has $c^+ \cdot d^+ \ll |\delta| \cdot \alpha^*$, the worst-case complexity is then in $O(|Q| + |\delta| \cdot \alpha^* + M \cdot c^+ \cdot d^*)$. A more precise analysis requires the definition of a notion of length associated with terms, opening the way for further optimizations, e.g. as described by Example 9 in the context of RNA folding.
### Example 9: RNA BP folding – Statistical sampling

In the context of RNA conformational sampling, Algorithm 2 specializes into the following backtrack, preceded by a computation of the partition function $Z_{[i,j]}$ for all $[i, j] \in [1, n]$, as described in Example 5.

**Function** SampleRNA2D ([i, j] ⊆ [1, n], w ∈ {A, C, G, U}ⁿ, Z):

if j - i ≤ θ then return ∅;

▷ Interval too short to support base pair

r ← UnifRand(Z_{[i,j]})

▷ Random r ∈ [0, Z_{[i,j]}[

r ← r - Z_{[i+1,j]};

if r < 0 then

| return SampleRNA2D ([i + 1, j], w)

for k ← i + θ + 1 to j do

if wᵢ, wₖ ∈ B then

r ← r - eβ × Z_{[i+1,k-1]} × Z_{[k+1,j]};

▷ Defining Z_{[j+1,j]} ≡ 1

if r < 0 then

return

{{i, k}} ∪ SampleRNA2D([i + 1, k − 1], w) ∪ SampleRNA2D([k + 1, j], w)

This function is a strict specialization of Algorithm 2, and randomly generates a secondary structure from the Boltzmann distribution when invoked on [1, n] with a sequence of length n.

As previously analyzed [148], its execution requires Θ(n²) operations in the worst-case scenario due to the possible Θ(n) iterations of the **for** loop, followed by a recursion on a subinterval only marginally smaller (n − 2). Its average-case behavior in the homopolymer model (all bases allowed to pair), or on a random RNA sequence of length n, was shown to be in Θ(n√n), a behavior that holds for a large class of combinatorial classes [71].

The time complexity can be dramatically lowered to O(n log n) through a minor modification of the algorithm. Indeed, it suffices to substitute a **Boustrophedon order**

$$i + \theta + 1 \sim j \sim i + \theta + 2 \sim j - 1 \sim i + \theta + 3 \sim \dots \sim \left\lfloor \frac{i + \theta + 1 + j}{2} \right\rfloor$$

instead of the (implicit) sequential order i + θ + 1 ∼ i + θ + 2 ∼ ... ∼ j within the **for** loop to instantly obtain a O(n log n) time complexity for the backtracks in the worst-case scenario.

### 2.3.3 Non redundant sampling (118, 133, 165)

As shown in Section ??, the level of redundancy can be overwhelming within a Boltzmann-Gibbs sample. One the one hand, redundancy is instrumental to the **consistency** of the estimator (2.16), since convergence towards the correct expectation requires that the frequency of any given element converges towards its probability. On the other hand, redundancy appears to be **non-informative**, and even **wasteful**, while sampling from a Boltzmann-Gibbs distribution, or any distribution known *a priori*. Indeed, the exact emission probability of a generated sample can be derived exactly from the partition function and, in principle, should not need to be estimated from the frequency. Moreover, sampling can be used to recover diverse dominant solutions, with no further involvement in statistical estimates, e.g. in RNA kinetics studies [110, 133] or in automated software
testing [189].

One possible way to avoid redundancy is to perform **sampling without replacement**, forbidding the generation at the m-th step of any previously generated element [s<sub>1</sub>, ..., s<sub>m-1</sub>]. The targeted emission probability P(s<sub>m</sub> | s<sub>1</sub>, ..., s<sub>m-1</sub>) is then

$$P(s_m | s_1, \dots, s_{m-1}) = \frac{P(s_m)}{1 - \sum_{i=1}^{m-1} P(s_i)} \quad (2.18)$$

where P(s) denotes the classic (redundant) Boltzmann-Gibbs probability. Note that successive generations are now dependent, so classic estimators such as the empirical mean of Equation (2.16) become biased and should not be used.

However, as shown by Rovetta et al [165], the expectation E(F(S)) of any given function F : Ω → ℝ in the Boltzmann distribution (or any known distribution) can still be estimated from a non-redundant sample [s<sub>1</sub>, ..., s<sub>m</sub>] by

$$\tilde{F}([s_1, \dots, s_m]) = \frac{1}{m} \sum_{i=1}^{m} F(s_i) \left( 1 + (m-i) \times P(s_i) - \sum_{j=1}^{i-1} P(s_j) \right).$$

This estimator is provably unbiased, consistent and always yields lower expected variance that the empirical mean, motivating the following algorithmic problem.

### Problem 8 (Non-redundant sampling):

**Input:** Unambiguous DP scheme Δ + scoring f, correct w.r.t. energy function E; β ∈ ℝ a constant; M a **number of samples**

**Output:** A **non-redundant sequence** [s<sub>1</sub>, ..., s<sub>M</sub>] ∈ Ω<sub>Δ</sub><sup>M</sup>, of elements:

$$P(s_m | s_1, \dots, s_{m-1}) = \frac{P(s_m)}{1 - \sum_{i=1}^{m-1} P(s_i)}$$

In collaboration with Andy Lorenz [118], we considered a restriction of the above problem to the language of a weighted context-free grammar. Its main idea is to reinterpret derivation as a **purely sequential process**, where derivations involved in the production of a term are transformed into a **sequence of derivations** through preorder traversal of the derivation tree. Below is an example of a derivation tree t and its linear representation ω<sub>t</sub>:

$$t = \frac{d_1}{d_2 \quad d_3 \quad d_6} \quad \Leftrightarrow \quad \omega_t = d_1 \cdot d_2 \cdot d_3 \cdot d_4 \cdot d_5 \cdot d_6$$

Note that the derivation tree t can be unambiguously recovered from a sequence ω<sub>t</sub> through a prefix evaluation, and reinterpreted as a term T<sub>ω<sub>t</sub></sub>.
<caption>**Algorithm 3:** Non-redundant sampling algorithm</caption>

**Input:** Unambiguous DP scheme $\Delta$; Number $M$ of samples.
**Output:** Set of $M$ random terms, distributed according to the non-redundant distribution of Equation (2.18).

**Compute $\mathcal{Z}_q, \forall q \in Q$**
*   ▷ Partition function algorithm (c.f. Problem 2)
*   $P_F \leftarrow (\perp, \emptyset)$; ▷ Tree initially restricted to root $\perp$
*   $\pi(\perp) \leftarrow 0$;
*   $\mathcal{F} \leftarrow \emptyset$ ▷ Terms generated over previous iterations $\rightarrow$ avoided

**for** $i \in [1, M]$ **do**
| $(\omega, Z_{\omega}, u_{\omega}) \leftarrow (q_{\omega}, \mathcal{Z}_{q_{\omega}}, \perp)$; ▷ Start from root
| $\mathcal{N} \leftarrow \{u_{\omega}\}$ ▷ Nodes traversed during generation
**while** $\omega \cap Q \neq \emptyset$ **do**
|   $\omega_1 \cdot q_L \cdot \omega_2 \leftarrow \omega$ s.t. $\omega_1 \in \delta^*$ and $q_L \in Q$ ▷ $q_L$ lefmost state in $\omega$
|   $r \leftarrow \text{UnifRand}(Z_{\omega} - \pi(u_{\omega}))$ ▷ Random $r \in [0, \mathcal{Z}_{\omega} - \overline{\mathcal{Z}}(\omega)]$
|   **foreach** $d := q_L \xrightarrow{\lambda} q_1 \cdots q_k \in \delta$ **do**
|   |   $\omega' \leftarrow \omega_1 \cdot d \cdot q_1 \cdots q_k \cdot \omega_2$; ▷ Simulate derivation $d$ in $\omega$
|   |   $u_{\omega'} \leftarrow \text{Child}(P_F, u_{\omega}, d)$; ▷ Creates $u_{\omega'}$ if absent $(\pi(u_{\omega'}) \leftarrow 0)$
|   |   $Z_{\omega'} \leftarrow Z_{\omega} \times (e^{-\beta \cdot E(\lambda)} \prod_{j=1}^{k} \mathcal{Z}_{q_j}) / \mathcal{Z}_{q_L}$; ▷ Update $Z$ per (2.19)
|   |   $r \leftarrow r - (Z_{\omega'} - \pi(u_{\omega'}))$;
|   |   **if** $r < 0$ **then**
|   |   |   ▷ Happens with probability $\frac{\mathcal{Z}_{\omega'} - \overline{\mathcal{Z}}(\omega')}{\mathcal{Z}_{\omega} - \overline{\mathcal{Z}}(\omega)}$
|   |   |   ▷ Move to selected child
|   |   $(\omega, Z_{\omega}, u_{\omega}) \leftarrow (\omega', Z_{\omega'}, u_{\omega'})$; $\mathcal{N} \leftarrow \mathcal{N} \cup \{u_{\omega'}\}$;
|   |   **break;**
|   $\mathcal{F} \leftarrow \mathcal{F} \cup \{T_{\omega}\}$; ▷ Interpretation of $\omega \in \delta^*$ as a term $T_{\omega}$
|   $\pi_{\omega} \leftarrow \prod_{\lambda \in \omega} e^{-\beta \cdot E(\lambda)}$;
**foreach** $u \in \mathcal{N}$ **do** $\pi(u) \leftarrow \pi(u) + \pi_{\omega}$; ▷ Update traversed node weights
**return** $\mathcal{F}$; ▷ Final non-redundant list of generated elements

performed during the generation $\mathcal{F}$. Any **node** $u_{\omega} \in P_F$ represents (implicitly) some immature word $\omega \in \mathcal{D}_F$. There exists a **directed edge** $(u_{\omega} \rightarrow u_{\omega'}) \in P_F$, **labeled by a derivation** $d(u_{\omega} \rightarrow u_{\omega'}) \in \delta$ if and only $\omega \xrightarrow{d} \omega'$ occurs in some sequence of $\mathcal{D}_F$. Finally, on each **node** $u \in \mathcal{D}_F$ we store a weight $\pi(u)$, updated after each **generation** such that the invariant $\pi(u_{\omega}) = \overline{\mathcal{Z}}(\omega)$ is maintained, allowing $\mathcal{O}(1)$ access to $\overline{\mathcal{Z}}(\omega)$ while computing probabilities (2.20).

We obtain Algorithm 3, which generalizes earlier context-free versions [118, 133, 147] and can be essentially found in Juraj Michalík's PhD thesis [132]. The complexity of the algorithm is the same, up to implementation constants, as the redundant stochastic backtrack.
**Remark 2.3.1 (Rejection-based non-redundant sampling):** In the **uniform distribution** ($\beta = 0$), a rejection-based sampling constitutes a reasonable alternative to produce M unique elements. Indeed, the search space is typically exponential $|\Omega| \approx \zeta^n$ for some $\zeta > 1$. A classic **coupon collector analysis** shows that the expected number of generations to get $|\Omega|$ distinct elements, is $|\Omega| \cdot H_{|\Omega|}$ with $H_m = \sum_{i=1}^m 1/i \in \Theta(\log n)$ the m-th harmonic number. Moreover, the expected number of generation is an increasing function of M, and is bounded by $M \log |\Omega| \in \Theta(M \cdot n)$. In other words, a repeated generation of elements until M distinct elements are produced has $\Theta(M \cdot n)$ expected time in a uniform distribution, i.e. non-redundancy sampling only induces a **linear overhead** in comparison with the redundant one.

In Boltzmann-Gibbs sampling, however, the overhead is expected to grow exponentially with n, so a **rejection-based strategy is exponentially less efficient than Algorithm 3**.

### 2.3.4 Adaptive sampling of constrained sequences (18, 90, 158)

As shown in Section 2.2.5, using explicit convolutions to compute the classified partition function quickly becomes prohibitively costly when (combinations of) expressive scoring functions are considered. Unfortunately, such a computation would be required to perform a **constrained sampling**, by adapting the stochastic backtrack of Section 2.3.2 to restrict the generation to objects having a value of interest.

As an alternative, in collaboration with Olivier Bodini [18], we introduced a **multidimensional Boltzmann sampling** method, also called **adaptive sampling** in the context of Bioinformatics applications [193] (due to the ubiquity of Boltzmann in Bioinformatics, leading to an overloaded nomenclature). This method is inspired both by the versatility of rejection methods in random generation (see Section 4.1), and by the typical concentration of distributions for additive parameters. This concentration is illustrated in 1D by Figure 2.1 for the number of base pairs in RNA, and in 2D by Figure 2.2 for two types of base pairs distinguishable in the basic Nussinov DP scheme.

**Problem 9 (Constrained sampling):**

**Input:** Unambiguous DP scheme $\Delta$ + scoring f, correct w.r.t. function E; $\beta \in \mathbb{R}$; **scoring functions** $F_1, \dots, F_k : \Lambda \to \mathbb{N}^+$ with $V := \max_{T \in \mathcal{T}_\Delta} F(T)$ associated with **objective values** $v_1^*, \dots, v_k^* \in \mathbb{N}^k$; and **number of samples** M

**Output:** A collection $C := [s_1, \dots, s_M]$ of random, Boltzmann-distributed w.r.t. E, elements of $\Omega_\Delta$ such that

$$F_i(s_j) = v_i^*, \forall s_j \in C, \forall i \in [1, k]. \qquad (2.21)$$

Rather than computing the costly convolution products described in Equation (2.12), we introduced a rejection based approach called **multidimensional Boltzmann sampling**,
<figure><img src="image_20.png" /><figcaption>Figure 2.2: **Impact of weights on the distribution of base pairs in RNA 2D folding.** We distinguish **maybe-stacked** base pairs (MS), i.e. involving both ends (i, j) of an investigated interval [i, j] (alt. such that j + 1 is base paired) from other base pairs, called **unstacked** (US). Those two contributions can be distinguished on the classic Nussinov-Jacobson DP scheme, such that a weighted distribution can be induced.</figcaption></figure>

whose principles are further described in Section 4.2.2. At its core, one uses a **weighting scheme** where numerical values ($\pi_1, \dots, \pi_k$) are associated to scoring functions, inducing a **weighted distribution**

$$P_{\pi_1, \dots, \pi_k}(T) = \frac{e^{-\beta \cdot E(T)} \prod_{i=1}^{k} \pi_i^{F_i(T)}}{\mathcal{Z}_{\pi_1, \dots, \pi_k}} \quad (2.22)$$

where $\mathcal{Z}_{\pi_1, \dots, \pi_k}$ is the weighted equivalent of the partition function.

Sampling within a weighted distribution can be done using Algorithm 2 based on a unifying scoring function $E^*$ that aggregates all the contributions of individual features, such that

$$E^*(T) := \sum_{\lambda \in T} E(\lambda) - \frac{\sum_{i=1}^{k} F_i(\lambda) \times \log \pi_i}{\beta} = E(T) - F_i(T) \times \frac{\sum_{i=1}^{k} \log \pi_i}{\beta}.$$

This generation is coupled with a **rejection** step, which filters out the objects that do not fulfill the constraints described in Equation (2.21).

A crucial observation is that, in many cases, the successive choices performed during a stochastic backtrack can be interpreted as a large collection of **independent events** (i.e. occurrences of constructors), associated with bounded values of which the scoring functions are just the sum. The **central limit theorem** thus applies, and induces distributions that are (multidimensional) normal (a.k.a. Gaussian) with relatively low variance/co-variance values. In particular, once the weights are **calibrated** such that

$$\mu_i := \mathbb{E}_{\pi_1, \dots, \pi_k} (F_i(S)) = v_i^*, \forall i \in [1, k],$$
<figure><img src="image_21.png" /><figcaption>Figure 2.3: Using Algorithm 4 to sample 10 000 Boltzmann-distributed secondary structures having exactly 15 base pairs. For the RNA w = (ACGU)²⁵ the expected number of base pairs is approximately 32.51 in the Boltzmann distribution, as shown in Figure 2.1 (T = 37°C, E(S) = -|S|). A binary search is used to determine a weight π<sub>BP</sub> such that the expected number of base pairs, empirically estimated at each iteration, matches the targeted one, accumulating across iterations.</figcaption></figure>

then one has

$$ \mathbb{P}(F_i(S) = v_i^*, \forall i \in [1, k]) = \frac{1}{\sqrt{(2\pi)^k \det \Sigma}} $$

where Σ is the co-variance matrix induces by scoring functions, such that Σ<sub>i,j</sub> := E((F<sub>i</sub>(S) - μ<sub>i</sub>)(F<sub>j</sub>(S) - μ<sub>j</sub>)). Assuming that k is constant, and denoting by V<sub>i</sub> the max value taken by F<sub>i</sub> over Ω<sub>Δ</sub>, one expects Θ(M√Π<sub>i</sub> V<sub>i</sub>) attempts in order to obtain M suitable element.

Our overall strategy to solve Problem 9 is summarized by (pseudo) Algorithm 4, and performs the weights optimization in parallel to the generation, accumulating suitable samples. Indeed, while following different overall distributions, samples generated across iterations of the weight optimization still follow the distribution (2.21) after the rejection step (see Section 4.1 and 4.2.1 for further discussions). The **weight optimization** can usually be tackled efficiently, since the expected number of occurrences of a given constructor is a strictly increasing function of the associated weight. The problem can generally be solved using tools in convex optimization [15], but can also be solved using **steepest descent** heuristics, or even a simple **binary search** in a mono-objective setting, are usually sufficient to solve the problem.

### Example 10: RNA BP folding – Probing the Boltzmann ensemble of higher energy

Let us illustrate the general principle, and typical complexity analysis, of multidimensional sampling in a simple setting. Here, we want to estimate statistical properties of 2D structures with exactly v<sub>BP</sub><sup>*</sup> base pairs, potentially far from the MFE, energy-wise, in a base-pair based energy model. First, we introduce a weight W<sub>BP</sub> ∈ ]0, +∞[ that controls the expected number of base-pairs. Its incorporation to the partition function computation (see Example 5) and sampling procedure (see Example 9) only requires the introduction of π<sub>BP</sub> within the contribution of valid base pairs, which becomes π<sub>BP</sub> × e<sup>β</sup> ×
**Algorithm 4:** High-level description of a multidimensional Boltzmann strategy for constrained sampling

**Input:** Num. samples $M$; Functions $F_1, \dots, F_k$ with targeted values $(v_i^*)_{i=1}^k$;

Estimate $\gamma$ for $\sqrt{(2\pi)^k \det \Sigma}$

**Output:** Collection $\{s_j\}_{j=1}^M$ of random, Boltzmann-distributed, elements such that $F_i(s) = v_i^*, \forall i \in [1, k]$

$C \leftarrow \emptyset$;

$\pi \equiv (\pi_1, \dots, \pi_k) \leftarrow (1, \dots, 1)$;

**while** $|C| < M$ **do**

$D \leftarrow$ Sample $M \times \gamma$ elements in the $\pi$-weighted distribution;

Update weight vectors towards $E_{\pi_1, \dots, \pi_k}(F_i(S)) = v_i^*, \forall i \in [1, k]$;

$C \leftarrow C \cup \{s \in D \mid F_i(s) = v_i^*, \forall i \in [1, k]\}$;

**return** $C$;

$$Z_{[i+1,k-1]} \times Z_{[k+1,j]}.$$

Then a simple binary search, explicit in the pseudocode below, is typically sufficient to guarantee a quick convergence towards suitable weights, as illustrated by Figure 2.3, and we get

$C \leftarrow \emptyset$;

($\pi_{BP}, \pi_{min}, \pi_{max}) \leftarrow (1, 0, +\infty)$;

**while** $|C| < M$ **do**

Compute $Z_{BP}$;

$D \leftarrow \{\text{SampleRNA2D}([1, n], w, Z_{BP})\}_{i=1}^{2M}$;

$\overbrace{\#bp} \rightarrow \frac{\sum_{S \in D} |S|}{|D|};$

if $\#bp < v_{BP}^*$ then

($\pi_{BP}, \pi_{min}, \pi_{max}) \leftarrow (\min(2\pi_{BP}, \frac{\pi_{BP}+\pi_{max}}{2}), \pi_{BP}, \pi_{max})$;

else

($\pi_{BP}, \pi_{min}, \pi_{max}) \leftarrow (\frac{\pi_{min}+\pi_{BP}}{2}, \pi_{min}, \pi_{BP})$;

$C \leftarrow C \cup \{S \in D \mid |S| = v_{BP}^*\}$;

**return** $C$;

▷ Using $\Theta(n^3)$ algorithm

▷ Estimate expected #BPs for $\pi_{BP}$

▷ Binary search update

## 2.4 Applications of ensemble dynamic programming

As can be seen in the previous section, the design of a suitable dynamic programming algorithm unlocks a wide array of potential inquiries regarding the properties of the ensemble of solutions, using algorithms that can be generically derived. As a consequence, the most crucial aspect of a DP algorithm, arguably the only one worthy of human attention, is the **design of a suitable dynamic programming scheme**, and the combinatorial decomposition onto which it is based. In particular, it is not only essential
<table><thead><tr><th>Conformation space</th><th>Base-pairs</th><th>Stacking-Pairs</th><th>Nearest-Neighbor</th></tr></thead><tbody><tr><td>Non-crossing</td><td>MFE<br>Approx.</td><td>P<br>Nussinov and Jacobson [143]<br>—</td><td>P<br>Leong et al. [100]<br>—</td><td>P<br>Zuker and Stiegler [213]<br>—</td></tr><tr><td>Planar</td><td>MFE<br>Approx.</td><td>???<br>2-approx.<br>±Leong et al. [100]</td><td>NP-Hard<br>Leong et al. [100]<br>2-approx.<br>Leong et al. [100]</td><td>NP-Hard<br>Leong et al. [100]<br>???</td></tr><tr><td>General</td><td>MFE<br>Approx.</td><td>P<br>Tabaska et al. [181]<br>—</td><td>NP-Hard<br>Lyngsø [124], Sheikh et al. [175]<br>(any* Δ model)<br>ε-approx. ∈ O(n<sup>4/3</sup>)<br>Lyngsø [124], Sheikh et al. [175]<br>1/5 (any Δ model)</td><td>NP-Hard<br>Lyngsø and Pedersen [123], Akutsu [3]<br>APX-Hard<br>Sheikh et al. [175]</td></tr></tbody></table>

<table><caption>Table 2.2: Computational complexity of RNA 2D prediction with pseudoknots in the energy minimization paradigm.</caption><thead><tr><th>Conformation space</th><th>Base-pairs</th><th>Stacking-Pairs</th><th>Nearest-Neighbor</th></tr></thead><tbody><tr><td>Non-crossing</td><td>MFE<br>Approx.</td><td>P<br>Nussinov and Jacobson [143]<br>—</td><td>P<br>Leong et al. [100]<br>—</td><td>P<br>Zuker and Stiegler [213]<br>—</td></tr><tr><td>Planar</td><td>MFE<br>Approx.</td><td>???<br>2-approx.<br>±Leong et al. [100]</td><td>NP-Hard<br>Leong et al. [100]<br>2-approx.<br>Leong et al. [100]</td><td>NP-Hard<br>Leong et al. [100]<br>???</td></tr><tr><td>General</td><td>MFE<br>Approx.</td><td>P<br>Tabaska et al. [181]<br>—</td><td>NP-Hard<br>Lyngsø [124], Sheikh et al. [175]<br>(any* Δ model)<br>ε-approx. ∈ O(n<sup>4/3</sup>)<br>Lyngsø [124], Sheikh et al. [175]<br>1/5 (any Δ model)</td><td>NP-Hard<br>Lyngsø and Pedersen [123], Akutsu [3]<br>APX-Hard<br>Sheikh et al. [175]</td></tr></tbody></table>

to ensemble analysis that the DP schemes are correct, but also that they are unambiguous and complete. In the following, I briefly mention some contributed schemes, addressing particular problems in RNA bioinformatics and comparative genomics.

## 2.4.1 (Recursive) simple type pseudoknots and kissing hairpins (150)

The absence of crossing base pairs in ensemble of conformations explored by classic algorithms represents a real limitation, and can be mainly attributed to:

1. The computational complexity of folding with pseudoknots in realistic models [3, 123], as summarized in Table 2.2, which we extended to large classes of energy models in collaboration with Saad Sheikh and Rolf Backofen [175];

2. The difficulty to imprint, at a combinatorial level, complex geometric constraints (e.g. chain closure) induced by base pairs, leading to ensemble of conformations that are highly dominated by unrealistic structures.

To work around such issues, a very rich body of literature [3, 19, 29, 53, 154, 156, 164, 183] has been dedicated to the design of dynamic programming algorithms capturing a subset of feasible conformations. However, many of those algorithms relied on an ambiguous DP scheme, preventing the natural exploration of ensembles of structures.

In collaboration with Cédric Saule [150], we revisited or designed popular DP schemes, visually described in Figure 2.4, for RNA folding with pseudoknots. In order to describe the generic derivation of ensemble algorithms (partition function, statistical sampling,
A – Unambiguous/complete Θ(n⁴) version of a DP scheme for **simple type pseudoknots** (3)

<figure><img src="image_22.png" /></figure>

B – Unambiguous/complete Θ(n⁵)/Θ(n⁴) version of a DP scheme for recursive kissing hairpins (183)

<figure><img src="image_23.png" /><figcaption>Figure 2.4: Unambiguous DP schemes for RNA folding in presence/absence of (recursive) pseudoknots [150].</figcaption></figure>

inside/outside...), we used an **oriented hypergraph** (a.k.a. directed monotonic hypergraph) formalism introduced by Finkelstein and Roytberg [66], also occurring as a generic representation of DP in the field of natural language processing [98].

In addition to designing new DP schemes for two popular families, we adopted an enumerative combinatorics perspective to help in the proof of their unambiguity/correctness. Indeed, many conformation spaces amenable to dynamic programming can be enumerated, and generating functions $PK(z) = \sum_{n \ge 0} p_n z^n$ with $p_n$ the number of (pseudokotted) conformations over n nucleotides, have been established for the most popular DP schemes [168]. On the other hand, a full relaxation of base pair compatibility constraints (a.k.a. homopolymer model), allows to associate a generating function
$S_Δ(z) = \sum_{n \ge 0} |\Omega_{Δ_n}|, z^n$ to a family of DP schemes $Δ_1, Δ_2, ...$ for instances of increasing length. In particular, the equality $PK(z) = S_Δ(z)$, in conjunction with a proven completeness of the DP scheme, is sufficient to conclude on its unambiguity. Conversely, the conjunction of $PK(z) = S_Δ(z)$ and a proven unambiguity implies the completeness, effectively halving the work required to prove, typically through tedious inductions, those two properties. We used this technique to verify the properties of two contributed schemes described in Figure 2.4.

## 2.4.2 Dual partition functions and evolutionary robustness of RNAs

Over several collaborative projects, we considered a dual version of the partition function, integrating over the space of compatible sequence/structure pairs [192, 193], or over the set of sequences associated with one [117, 157, 158] or several [88, 90] given structures.

In a first collaboration with Jérôme Waldispühl [192, 193], we leveraged a dynamic programming scheme introduced by the RNAMutants [191] method, to investigate the interplay between sequence and structure at prescribed GC-content. Its DP scheme explores all pairs $(w, S)$ of compatible sequence and structure, classified by Hamming distance $δ(w, w^*)$ to a reference wild type sequence $w^*$ of length $n$, allowing for a $Θ(n^5)/Θ(n^3)$ computation of the partition function of k-mutants

$$Z_{[k]} := \sum_{w \in \Sigma^n \text{ s. t. } S \text{ comp. with } w} \sum_{\delta(w, w^*) = k} e^{-\beta \cdot E(w, S)}$$

A restricted version of the partition function, only considering the subinterval $w^*_{i,j}$ of sequence $w^*$, can then be computed through

$$Z_{i,j} = \sum_{[k] \in \Sigma} Z_{[k-\delta_{w^*_i, b}]} + \sum_{b, b' \in \Sigma^2} \sum_{l=i+\theta+1}^{j} \sum_{k'=0}^{k-\delta_{w^*_i, w^*_{l', b'}}} e^{-\frac{\beta \cdot E_{b, b'}}{RT}} \cdot Z_{[k]} \cdot Z_{[k-k'-\delta_{w^*_i, w^*_{l', b'}}]} \quad (2.23)$$

in the case of a base-pair based energy model (although the implementation supports the Turner nearest-neighbor energy model [185]).

The same partition function was used to perform statistical sampling, revealing a major caveat of the model: due to the overwhelming impact of the GC-content on the partition function of a given sequence, the Boltzmann distribution was observed to disproportionately focus on sequences having extreme GC-contents, potentially impacting the accuracy of conclusions drawn while studying the evolutionary robustness of concrete sequences. To circumvent this issue, we introduced an adaptive sampling algorithm, an instance of multidimensional Boltzmann sampling [18], to direct the sampling towards areas of predetermined GC-content, allowing us to study the interplay between evolutionary robustness and thermodynamic stability [192, 193].
<figure><img src="image_24.png" /><figcaption>Case 1: First position is unpaired.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Case 2: Extremities are paired, nested within a consecutive base-pair, forming a stacking base-pair.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Case 3: First position in paired to some position, but not involved in a stacking pair.</figcaption></figure>

Figure 2.5: Classified DP scheme for the mutational landscape with respect to a single dominant conformation [157].

Another instance of a DP scheme generating sequences considers the space of sequences compatible with a single structure, in an attempt to use structure conservation to detect **sequencing errors in metagenomic data sets** [157, 159]. In collaboration with Jérôme Waldispühl and Vladimir Reinharz, we considered an energy model of intermediate complexity, mainly including stacking pairs, and computed the **dual partition function** associated with a reference structure $S^*$ and sequence $w^*$:

$$Z_{S^*,[k]} = \sum_{\substack{w \text{ comp. with } S^* \\ \delta_{w,w^*} = k}} e^{-\beta \cdot E(w,S^*)}$$

using a modified version of (2.23), showcased in Figure 2.5, which runs in $\Theta(n \cdot M^2)$ time and $\Theta(n \cdot M)$ memory, with M the maximum number of suspected errors. In combination with an **outside contribution**, computed in a similar fashion, we were able to compute individual probabilities of assignments for each position. We also included contributions in our scoring scheme for the **isostericity** of base pair substitutions [115, 178] revealing a good capacity to discriminate, and to some extent correct, erroneous positions.

Dual partition functions are also essential for the **positive design of RNAs**, which we addressed as a random generation task over a series of collaborative works [88, 90, 117, 160, 209, 210]. In this context, the control of the maximum number of mutations is lifted, allowing a $\Theta(n)$ time/space computation in the case of a target secondary structure. Recently [88, 90], we tackled the problem in its most general setting, adopting
a declarative approach where constraints are associated to a subset of positions within an RNA, either restricting their potential joint content or inducing specific free-energy contributions. The entire network of dependencies induced by the constraints can then be abstracted as an hypergraph, whose **tree decompositions** can be used to derive systems of DP equations, optimal in some sense, without any direct human intervention! We refer the reader to Chapter 3 for further context and motivation.

### 2.4.3 Unambiguous tree alignments (32, 33)

Generalizing the notion of sequence alignment, a **tree alignment** for two trees $T_1$ and $T_2$, is a set of pairwise correspondences, or **matches**, between the nodes of $T_1$ and $T_2$, consistent with the ancestry relationships in $T_1$ and $T_2$. Nodes without a partner in $T_1$ (resp. $T_2$) are called insertions (resp. deletions) Equivalently, a tree alignment can be represented by a **supertree**, a tree $T$ whose nodes are pairs of nodes $(u_1, u_2) \in (T_1 \times T_2) \cup (T_1 \times \{-\}) \cup (\{-\} \times T_2)$ such that $T_1$ (resp. $T_2$) can be recovered from $T$ by considering the first component $u_1$ (resp. second component $u_2$) in each pair, and finally contracting insertions/deletions nodes $-$. Tree alignments are natural in the context of RNA bioinformatics, where they represent a reasonable approach to measure the similarity between two secondary structures [17, 171].

Algorithmically, a classic problem is to find the optimal alignment between two given trees $T_1$ and $T_2$, e.g. the supertree that minimizes the accumulated log odds of matches and indels. It is solved by Jiang, Wang and Zhang [103] (JWZ) based on the DP scheme:

$$ \text{Align} \left( \begin{array}{ccc} \text{Redundancy} & & \\ & \ddots & \\ & & \text{Redundancy} \\ \text{Redundancy} & & \\ & & \ddots \\ & & & \text{Redundancy} \\ & & & & \ddots \\ & & & & & \text{Redundancy} \\ & & & & & & \text{Redundancy} \end{array} \right) = \min \left\{ \begin{array}{l} \text{Align} \left( \begin{array}{ccc} \text{Redundancy} & & \\ & \ddots & \\ & & \text{Redundancy} \\ \text{Redundancy} & & \\ & & \ddots \\ & & & \text{Redundancy} \\ & & & & \ddots \\ & & & & & \text{Redundancy} \end{array} \right) + \text{Del}(\bullet) \\ \text{Align} \left( \begin{array}{ccc} \text{Redundancy} & & \\ & \ddots & \\ & & \text{Redundancy} \\ \text{Redundancy} & & \\ & & \ddots \\ & & & \text{Redundancy} \\ & & & & \ddots \\ & & & & & \text{Redundancy} \end{array} \right) + \text{Ins}(\bullet) \\ \text{Align} \left( \begin{array}{ccc} \text{Redundancy} & & \\ & \ddots & \\ & & \text{Redundancy} \\ \text{Redundancy} & & \\ & & \ddots \\ & & & \text{Redundancy} \\ & & & & \ddots \\ & & & & & \text{Redundancy} \end{array} \right) + \text{Subst}(\bullet, \bullet) \end{array} \right\} \\ \text{Align} \left( \begin{array}{ccc} \text{Redundancy} & & \\ & \ddots & \\ & & \text{Redundancy} \\ \text{Redundancy} & & \\ & & \ddots \\ & & & \text{Redundancy} \\ & & & & \ddots \\ & & & & & \text{Redundancy} \end{array} \right) = \min \left\{ \begin{array}{l} \text{min} \left\{ \text{Align} \left( \begin{array}{ccc} \text{Redundancy} & & \\ & \ddots & \\ & & \text{Redundancy} \\ \text{Redundancy} & & \\ & & \ddots \\ & & & \text{Redundancy} \\ & & & & \ddots \\ & & & & & \text{Redundancy} \end{array} \right) + \text{Align} \left( \begin{array}{ccc} \text{Redundancy} & & \\ & \ddots & \\ & & \text{Redundancy} \\ \text{Redundancy} & & \\ & & \ddots \\ & & & \text{Redundancy} \\ & & & & \ddots \\ & & & & & \text{Redundancy} \end{array} \right) + \text{Del}(\bullet) \right\} \\ \text{min} \left\{ \text{Align} \left( \begin{array}{ccc} \text{Redundancy} & & \\ & \ddots & \\ & & \text{Redundancy} \\ \text{Redundancy} & & \\ & & \ddots \\ & & & \text{Redundancy} \\ & & & & \ddots \\ & & & & & \text{Redundancy} \end{array} \right) + \text{Align} \left( \begin{array}{ccc} \text{Redundancy} & & \\ & \ddots & \\ & & \text{Redundancy} \\ \text{Redundancy} & & \\ & & \ddots \\ & & & \text{Redundancy} \\ & & & & \ddots \\ & & & & & \text{Redundancy} \end{array} \right) + \text{Ins}(\bullet) \right\} \\ \text{Align} \left( \begin{array}{ccc} \text{Redundancy} & & \\ & \ddots & \\ & & \text{Redundancy} \\ \text{Redundancy} & & \\ & & \ddots \\ & & & \text{Redundancy} \\ & & & & \ddots \\ & & & & & \text{Redundancy} \end{array} \right) + \text{Align} \left( \begin{array}{ccc} \text{Redundancy} & & \\ & \ddots & \\ & & \text{Redundancy} \\ \text{Redundancy} & & \\ & & \ddots \\ & & & \text{Redundancy} \\ & & & & \ddots \\ & & & & & \text{Redundancy} \end{array} \right) \end{array} \right\} $$

which can be computed in time $\Theta(n_1 \cdot n_2 \cdot \max(n_1, n_2)^2)$. Moreover, as shown by Herrbach et al. [92], the expected time complexity of the JWZ algorithm is actually in $\Theta(n_1 \cdot n_2)$ for uniform random trees, the quartic complexity being circumscribed to trees with large degrees, thus having extremely low probability.

However, this decomposition scheme features a **high level of redundancy**, i.e. an alignment of $T_1$ and $T_2$ has, on average, an exponential number (on $|T_1| + |T_2|$) of associated supertrees generated by the above DP scheme. Worse, some alignments are associated
<figure><img src="image_27.png" /><figcaption>Figure 2.6: Unambiguous DP scheme for tree alignments.</figcaption></figure>

by a unique supertree, so the algorithm cannot be easily adapted into a tool for ensemble
analysis (e.g. to associate a notion of support for a given alignment).

With Cédric Chauve and Julien Courtiel [32, 33], we designed a **correct, unambiguous and complete DP scheme for tree alignments**, summarized in Figure 2.6. While highly technical in its design and proof, it preserves both the worst-case $\Theta(n_1 \cdot n_2 \cdot \max(n_1, n_2)^2)$ and average-case $\Theta(n_1 \cdot n_2)$ complexities of the JWZ algorithm. In fact, using analytic combinatorics, we could show that the average-case complexity of both algorithms, expressed as a function of the overall size $n := n_1 + n_2$ of trees, only grows like $\Theta(n\sqrt{n})$. This counter-intuitive complexity, remarkably much lower than the $\Theta(n^2)$ average-case complexity of sequence, is however only an artifact induced by the bimodal partition of length across a pair of trees. Indeed, upon generating a pair of trees ($T_1, T_2$) uniformly at random such that $|T_1| + |T_2| = n$, one typically gets a tree of size $\Theta(\sqrt{n})$ and a tree of size $\Theta(n)$, rather than two sequences of length $n/2$ in the case of sequences. One should therefore exercise caution when performing an average case analysis on pairs of objects, lest the final result may be more representative of the length distribution than of the underlying phenomenon...

The DP scheme could also be interpreted as a grammar generating a subset of supertrees,
bijectively associated with tree alignments, and we used analytic combinatorics techniques to derive asymptotic properties of alignments (e.g. average # matches, number...). The same could be done for the (redundant) family of supertrees produced by the JWZ DP scheme, and taking the ratio gave the expected number of supertrees, generated by JWZ for each alignment, confirming the exponential level of ambiguity.
Chapter 3

RNA design

The good performances of RNA structure prediction, especially at the secondary structure
level, suggest harnessing their underlying models and methods to perform a **rational**
**design of RNAs**, broadly defined here as the construction of novel nucleotide sequences
achieving one (or several) predefined function(s). In the context of my work, a strong
emphasis is put on the development of **computational methods for an automated RNA**
**design**.

The originality of my contributions, spanning several distinct collaborations, lies in the adoption of a **random generation perspective** over the positive and, to some extent, negative design of RNA. Additionally, I have attempted to revisit negative design in a minimal setting, amenable to exact (if partial) combinatorial analysis, uncovering a promising, yet unexploited, notion of approximate design, that could lead to future developments.

**Outline.** I first describe in Section 3.1 the motivation, and various use-cases for RNA design. Section 3.2 describes partial exact solutions for a combinatorial version of the NP-hard inverse folding. Finally, Section 3.3 presents an original contribution to positive design based on a constrained random generation, surprisingly providing a satisfactory initialization steps for negative design.

The following summarizes, and attempts to unify, contributions described within the
following list of articles published in journals and/or presented at conferences.

Associated contributions

A. Levin, M. Lis, Y. Ponty, C. W. O'Donnell, S. Devadas, B. Berger, and J. Waldispühl. A global sampling approach to designing and reengineering RNA secondary structures. Nucleic Acids Research, 40(20):10041–10052, Nov. 2012
Y. Zhou, **Y. Ponty**, S. Vialette, J. Waldispühl, Y. Zhang, and A. Denise. Flexible RNA design under structure and sequence constraints using formal languages. In *ACM-BCB 2013*, Bethesda, Washington DC, United States, Sept. 2013

V. Reinharz, **Y. Ponty**, and J. Waldispühl. A weighted sampling algorithm for the design of RNA sequences with targeted secondary structure and nucleotide distribution. *Bioinformatics*, 29(13):i308–15, July 2013

<figure><img src="image_28.png" /></figure>

V. Reinharz, **Y. Ponty**, and J. Waldispühl. A weighted sampling algorithm for the design of RNA sequences with targeted secondary structure and nucleotides distribution. In *ISMB/ECCB 2013*, Berlin, Germany, July 2013

Y. Zhang, **Y. Ponty**, M. Blanchette, E. Lecuyer, and J. Waldispühl. SPARCS: a web server to analyze (un)structured regions in coding RNA sequences. *Nucleic Acids Research*, 41(Web Server issue):W480–5, July 2013

M. Drory Retwitzer, V. Reinharz, **Y. Ponty**, J. Waldispühl, and D. Barash. incaRNAfbinv : a web server for the fragment-based design of RNA sequences. *Nucleic Acids Research*, 44(W1):W308–W314, 2016

J. Hales, A. Héliou, J. Manuch, **Y. Ponty**, and L. Stacho. Combinatorial RNA Design: Designability and Structure-Approximating Algorithm in Watson-Crick and Nussinov-Jacobson Energy Models. *Algorithmica*, 79(3):835–856, 2017

<figure><img src="image_29.png" /></figure>

J. Hales, J. Manuch, **Y. Ponty**, and L. Stacho. Combinatorial RNA Design: Designability and Structure-Approximating Algorithm. In *CPM 2015*, Ischia Island, Italy, June 2015

A. Churkin, M. D. Retwitzer, V. Reinharz, **Y. Ponty**, J. Waldispühl, and D. Barash. Design of RNAs: comparing programs for inverse RNA folding. *Briefings in Bioinformatics*, 19(2):350–358, Jan. 2018

S. Hammer, W. Wang, S. Will, and **Y. Ponty**. Fixed-parameter tractable sampling for RNA design with multiple target structures. *BMC Bioinformatics*, 20(1):209, Dec. 2019

<figure><img src="image_30.png" /></figure>

S. Hammer, **Y. Ponty**, W. Wang, and S. Will. Fixed-Parameter Tractable Sampling for RNA Design with Multiple Target Structures. In *RECOMB 2018*, Paris, France, 2018

H.-T. Yao, C. Chauve, M. Regnier, and **Y. Ponty**. Exponentially few RNA structures are designable. In *ACM-BCB 2019*, pages 289–298, Niagara-Falls, United States, 2019. ACM Press
## 3.1 Why do we design RNAs?

At a fundamental level, designing RNAs represents the **ultimate stress-test** for our understanding of how RNA folds and acts on its environment. In this setting, one designs RNA sequences expected to fold into a predefined structure with respect to a folding prediction model. Synthesizing the resulting sequences, and using experimental methods to verify the actual adoption of the desired structure, one either validates the model or reveals some of its flaws, fueling and prioritizing further developments. Indeed, misfolding designed RNAs reveal gaps in our energy models and descriptors of the conformation spaces used by predictive algorithms. A similar strategy can be more generally used to test functional hypotheses involving the structure of RNA.

Molecular design also represents one of the primitives of **synthetic biology**, and RNAs have been used as biosensors [65], regulators, nano-material. . . Some naturally-occurring RNAs are sufficiently stable to fold in a modular fashion, enabling a *copy/paste approach* that simply combines existing RNAs. However, such a strategy is hindered, in an *in vivo* context, by the competition of artificial and endogenous RNAs, and by the intrinsic difficulty to produce orthogonal constructs over a limited number of available architectures. A rational design, coupled with an experimental filtering phase, is thus likely to represent the method of choice for future endeavors in RNA-based synthetic biology.

At a (primarily) sequence-based level, design is essential for future developments of **RNA-based therapeutics**. For instance, the recent discovery of viable treatments based on RNA interference is fueled by an understanding of how small RNAs can interfere with selected messenger RNAs to activate or inhibit them. In this context, an optimization of the nucleotide sequence of small interfering RNAs, akin to a design task, is crucial to ensure its efficacy and selectivity, mitigating the risk of side-effects for the drugs. Similarly, the specificity of genetic contents targeted by CRISPR-based editing, an thus its limited toxicity, is ensured by a redesign of guide RNAs. More generally, the adoption of a stable structure is very often a prerequisite for the interaction of RNAs with selected proteins [45], and are therefore crucial for the functionality of designed RNAs in a cellular context.

Design also helps in the **search for homologous RNAs**. Indeed, in multiple RNA families, selective pressure mostly applies at the structure level, hindering the discovery of new occurrences of a given RNA genes within the same organism (paralogs) or across available genomes (orthologs). If a structural model is known, and if the number of identified homologs is limited, a natural approach is to enrich homologs with sequences designed as to fold into the shared structure, in order to cover a larger proportion of the (neutral) sequence space.
### 3.1.1 The potential of RNA design for statistical evolutionary studies

I believe that RNA design should be systematically used within families of structured RNAs to assess the **significance of observed phenomena**. Indeed, a popular strategy consists in analyzing multiple sequence alignments of RNAs, and interpret quantitative observations in the light of an accepted pressures. To that purpose, one uses statistical tests to refute the hypothesis that the observation is merely induced by established properties. For instance, in the context of messenger RNAs, codon bias and the flanking of ORFs by start/stop codons should be accounted for when drawing consequences from the distribution of k-mers.

However, in the context of structural RNAs, most analyses only indirectly capture the overall stability (free-energy) of transcripts through a **preservation of the dinucleotide content**, based on the observation that sequences sharing a dinucleotide composition tend to have similar free-energy in their most stable 2D conformation. However, nothing prevents alternative dinucleotides compositions from achieving the same expected free-energy, so the distribution considered by classic tests is arguably *too narrow* to convincingly refute the hypothesis that an observed phenomenon can be imputed to the overall MFE stability (rather than be indicative of a new, yet to be explained, selective pressure).

Moreover, most families of functional RNAs not only share a comparable stability, but also share the adoption of a common functional structure. In this context, the significance of a property (e.g. overrepresentation of a motif, hypothesized as involved in a molecular recognition mechanism) should be assessed within the set of sequences that are believed fold into the common structure (+ additional established properties). In other words, one should be able to **perform a random generation of RNA designs** in order to suitably assess a notion of significance.

### 3.1.2 The different flavors of RNA design

RNA design targets a wide diversity of biological functions and, as such, encompasses of wide array of computational tasks. Two dominant paradigms dominate current computational approaches.

#### 3.1.2.1 Positive design

**Positive design** can be loosely defined as focusing on to **propensity** of produced RNAs to achieve a certain function. In a structural context, positive design usually involves
generating one or several sequences having good affinity (i.e. stability ≈ minimum free-energy) towards a targeted structure. A more functional perspective will also consider the presence/absence of sequences motifs, or the affinity towards interacting with another molecule, as assessed by an energy model (e.g. RNA/DNA or RNA/RNA interactions) or a probabilistic model (e.g. protein/RNA interactions). This may include compositional biases, such as the GC-content.

### 3.1.2.2 Negative design

**Negative design** tasks loosely focus on the **specificity** of produced sequences, attempting to avoid unwanted functions, interactions, structures... From a computational perspective, such tasks are usually harder than the one induced by positive design. Indeed, they typically require solving an **inverse combinatorial problem**, generally formulated as follows.

**Problem 10 (Inverse combinatorial problem):**

**Input:** Expected output $y^*$ and objective function $f$.

**Output:** Instance $x$ such that

$$ \{y | f(y; x) = \max_{y'} f(y'; x)\} = \{y^*\}. $$

The difficulty of those problems stems from the potential difficulty of optimizing the objective function, so that the natural decision version of the problem may not even belong in NP. Moreover, even when the optimization (coupled with the production of a certificate for the unicity of the optimal solution) can be performed in polynomial time, the inverse problem can still be NP-hard [20].

In the context of structural RNA design, negative design is usually referred to as the **inverse folding problem**, and historically consists in producing a nucleotide sequence, predicted to (uniquely) fold into a targeted structure with respect to the Minimum Free-Energy (MFE) criterion. Variants of the inverse folding problem consider the minimization of various notions of **defects**, including the: **MFE defect**, the free-energy difference between the target and MFE structure (or first suboptimal structure if the target and MFE structures coincide); **probability defect**, the Boltzmann probability for the sequence to adopt an alternative structure; or the **ensemble defect**, the expected base-pair distance between the target and a random structure in the Boltzmann-Gibbs distribution.
### 3.1.2.3 Computational approaches

Due to the robust NP-hardness of inverse folding [20, 91, 172], the dominant paradigm, shared by multiple computational approaches [6, 8, 27, 89, 94, 115, 208], combine an initial generation of **seed sequences**, akin to a (sometimes very rudimentary) positive design, with further rounds of **optimization** using a **local search**. Inverse folding has also been tackled using a variety of metaheuristics, including genetic algorithms [63, 64, 125], ant-colony optimization [107]....

Finally, the imperfect nature of the objectives of RNA design has left the door open for machine learning, recently attracting the predatory attention of neural networks [166] and/or crowdsourcing [109], with modest performances and disputable rationales.

## 3.2 Exact combinatorial negative design (85, 86)

The **computational hardness of inverse folding**, which had essentially eluded a whole community for almost three decades, was finally rigorously established in 2018 by Bonnet et al. [20].

It was preceded by an effort worthy of a mention by Schnall-Levin et al. [172], who generalized the RNA design problem into the **inverse Viterbi problem**, consisting of figuring out a sequence w, in a probabilistic grammar such as a e.g. Hidden Markov Model (HMMs) or a Stochastic Context-Free Grammar (SCFG), such that a predefined sequence of states is the most likely parsing for w. The problem was found to be NP-hard, even for HMMs, leading the authors to argue on an implied hardness of RNA design, since RNA folding can be viewed as a special case of SCFG parsing. However, while the hardness of the inverse Viterbi problem for SCFGs clearly follows from their result, it did not imply the hardness of the problem for a given input grammar, as is the case in RNA design.

Over the course of a collaboration with Jozef Hales, Alice Héliou, Jan Manuch and Ladislav Stacho [85, 86], we investigated a minimal installment of negative design, aiming to provide partial characterization for the space of designable structures. Namely, we consider the **Watson-Crick energy model**, initially restricting admissible base pairs to {A, U} and {G, C}. In this restricted context, inverse folding greatly simplifies and we get the following formulation.
**Problem 11 (Watson-Crick combinatorial design):**

**Input:** Secondary Structure S* over n nucleotides

**Output:** Sequence w such that

$$ \forall S \in S_n \text{ compatible with } w, S \neq S^* : |S| < |S^*|. $$

Bonnet et al. [20] showed that a slightly restricted version of the problem, forcing the assignment of predefined nucleotides to positions in w, is NP-hard. We thus focused our efforts on characterizing of designable subclasses of structures.

### 3.2.1 Basic results

Our results require a reminder of the representation for the secondary structure, illustrated in Figure 3.1. Such a representation abstracts a structure as a tree, rooted in a virtual node, with other nodes either being **base-pairs nodes** (labeled by interacting positions) or **unpaired nodes** positions (labeled by a single position). We define the **inner degree** of a node as its number of base pair neighbors. For instance, the node 4-37 in Figure 3.1 has inner degree 3.

**Theorem 1 (Basic designable classes):** Any structure S, represented as a tree, can designed in linear time as soon as (at least one) of the following conditions is met:

1. Each node in S has inner degree at most two;

2. S is saturated (no unpaired base) + any node has inner degree at most four;

Moreover, we established that the set of designable structures is **closed under k-stuttering**. In other words, if a structure S, represented as a dot-parenthesis word $S_1.S_2...S_n$, admits a known design $w := w_1.w_2...w_n$, then

$$ w_1^k.w_2^k...w_n^k \text{ is a design for } S_1^k.S_2^k...S_n^k, \forall k > 0. $$

The stuttering operation can also be described at the tree level, and essentially (up to renumbering of nodes) amounts to: i) duplicating each unpaired node k times; and ii) replacing any base pair node with a sequence of base pair nodes of length k.

We also identified two undesignable motifs, a.k.a. **local obstructions**, whose presence is sufficient to forbid any design for the structure. Namely, it is probably impossible to design any structure containing a multiloop with degree $\ge 5$ (m_5 **motif**), or with degree $\ge 3$ with $\ge 1$ unpaired base(s) (m_3, **motif**). Such motifs are universal to any
<figure><img src="image_31.png" /><figcaption>Figure 3.1: **Tree representation and separated coloring for a designable secondary structure.** The secondary structure admits a separated black/white/gray coloring of its base pairs, i.e. levels (red callouts) of gray base pairs (1 and 0) and unpaired base/leaves (2 and 4) do not overlap. The structure can thus be designed, e.g. through the sequence below.</figcaption></figure>

version of design, as they are essentially a consequence of the local nature of classic energy models, defects definition and DP-schemes used for structure prediction. We exploited this property in subsequent work focusing on the enumeration of designable structures Yao et al. [206], further described in Section 5.2.4.

### 3.2.2 Design as a tree coloring problem

Finally, we introduced a sufficient condition for designability, based on a special coloring of the tree representation. Specifically, we considered a **3-coloring** $\chi$, associating either black (G-C), white (C-G) or gray (A-U or U-A) to the base-pair nodes of a structure S. Unpaired leaves are not colored, and are typically assigned a U. We define $\chi$ as **proper** such that:

* Any node has at most one black node, one white node and two gray nodes;

* White nodes have no black children;

* Black nodes have no white children;

* Gray nodes have at most one gray child.

Any failure to obey those rules implies the existence of **local rewirings of base pairs**, leading to alternative configurations, and the sequence implied by the coloring is thus not a design. However, a suitable design also requires the global absence of competing structures, so we need a second property to conclude on the validity of the induced design.
We introduce the notion of **level** level_χ(v) of a node v within a coloring χ, defined as the number of black nodes, minus the number of white nodes, encountered while going from the node to the root. It can also be defined recursively as

$$ \text{level}_{\chi}(v) = \begin{cases} 0 & \text{if } v = \text{Root} \\ \text{level}_{\chi}(p) + 1 & \text{if } \chi(p_{v}) = \text{black} \\ \text{level}_{\chi}(p) - 1 & \text{if } \chi(p_{v}) = \text{white} \\ \text{level}_{\chi}(p) & \text{if } \chi(p_{v}) = \text{gray.} \end{cases} $$

where p_v denotes the parent of v. Figure 3.1 illustrates the levels associated with a proper coloring.

**Definition 3.2.1 (Separated 3-coloring):** Given a secondary structure S, a proper 3-coloring χ is **separated** if and only if

$$ \{\text{level}(v) | v \text{ base paired in } S\} \cap \{\text{level}(v) | v \text{ unpaired in } S\} = \emptyset. $$

This property is actually sufficient to rule out the absence of an alternative structure having as many (or more) base pairs. The underlying intuition is that, in any competing structure, some U from an unpaired node must be paired to an A within a gray node. In a separated coloring, those interacting nucleotides are at different levels, so the alternative A = U base pairs provably delimits two regions, each having different numbers of Gs and Cs. It follows that, since pseudoknots are forbidden, then some G and C must remain unpaired in the final structure. In other words, the competing structure has less base pairs than the target.

**Theorem 2 (Designability of separable structures):** Any secondary structure S that admits a proper and separated 3-coloring χ can be designed. Moreover, if χ is known, then a design can be produced in linear time.

Figure 3.1 illustrates the concept of separated coloring, and some associated sequence. It is worth noting that the existence of a separated coloring does not represent a necessary condition for design. Moreover, finding a proper/separated coloring remains an open algorithmic problem, so this result may appear, at first sight, devoid of any practical consequence. Nonetheless, it is at the root of the notion of structure-approximating design, developed over the next section.
<figure><img src="image_32.png" /><figcaption>**Figure 3.2: Execution of the structure-approximating design algorithm.** An admissible target structure (1) is represented as a tree, and decorated into a proper 3-coloring by a greedy algorithm (2). A separated proper coloring (3) is obtained by adding selected nodes (dotted) to prevent direct adjacency and assign different parities for gray and unpaired nodes. The resulting coloring guarantees a linear-time design for a structure (4) that only differs by a few base pairs (purple) from the target.</figcaption></figure>

### 3.2.3 Structure-approximating design

While it remains unclear how a separated coloring can be efficiently found, it is easy to adapt a proper, yet non-separated, coloring into a separated one through the insertion of a limited number of nodes. Indeed, one can insert a small number of additional base paired nodes in order to separate gray paired and unpaired nodes, i.e. assigning them different levels to ensure the separated nature of the resulting coloring. This strategy works for any tree that avoids the $m_5$ and $m_{30}$ motifs, noted that the presence of those motifs remains unchanged by the addition of base pairs.

A first idea, not necessarily optimally parsimonious, consists in **offsetting gray nodes to even levels**, and **unpaired nodes to odd levels**. This can be done, for any structure S avoiding $m_5$ and $m_{30}$, in three steps:

1. Produce an initial **greedy proper coloring** of the tree:

* Assign to the ($\ge$ 4) children of the root, colors black→white→gray→gray in this order, and descend recursively into its children;

* On any internal nodes, colored in black (resp. white and gray), assign colors b→g→g (resp. w→g→g and b→w→g) to its children;

2. Correct the coloring into a separated one:

* Insert a black paired node to prevent any direct adjacency between a gray paired node and an unpaired leaf. Note that such cases only arise for gray nodes associated with **isolated base pairs**;

* Traversing the nodes of the tree in preorder, duplicate the parent of any gray node at odd level, or unpaired node at even level (the parent of such nodes is provably either black or white).
After executing this algorithm on a structure S, one obtains a tree/coloring pair (S', χ') where χ' is separated. Moreover, at most two base pairs are added to each of the helices of S (brought down to one base pairs in the absence of isolated base pairs).

**Theorem 3 (Approximated design):** Given a structure S (resp. with no isolated base pair) avoiding m₅ and m₃₀, then a designable structure S' can be obtained by adding ≤ 2 (resp. ≤ 1) base pair(s) to each helix in S. Moreover, S' can be built, and designed for, in linear time.

In practice, the number of added base-pairs usually remains modest, as shown in Figure 3.2. Future extensions of this work include the design of algorithms for adding a minimal collection of base pairs to make the structure/coloring separated.

More importantly, we believe that the **separated** property should find a natural equivalent in more sophisticated energy models, including as the Turner energy model [185]. A first extension of our main results to virtually any weighted base-pair model substantiates this hope.

**Proposition 4 (Extension to weighted models):** Theorems 1, 2 and 3 hold in any energy model that is defined additively on base pairs, each having individual free-energy contributions:

$$E(G, C) = E(C, G) \rightarrow \alpha, E(A, U) = E(U, A) \rightarrow \beta \text{ and } E(G, U) = E(U, G) \rightarrow \gamma$$

such that $\alpha, \beta < \gamma$.

## 3.3 Stochastic positive design under constraints

In the context of positive design, computational objectives are typically sufficiently simple to be amenable to a random generation. Such an approach is preferable to a **one-shot resolution** of the constraints imposed by the model of functions, for a variety of reasons:

* Positive design is typically used to generate an initial sequence, the **seed sequence** of existing methods based on meta-heuristics. Generating the initial at random thus helps meta-heuristics (local search, or genetic algorithms) overcome barriers in the objective function;

* Although currently underexploited, statistics derived from a set of random designs provide information on a null model for RNA evolution, allowing to assess the significance of observed phenomenon [151, 209];
<figure><img src="image_33.png" /><figcaption>1) Boltzmann probability vs Free-energy of target</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>2) GC-content vs BP distance between MFE and target</figcaption></figure>

Figure 3.3: **Comparison of uniform and Boltzmann-Gibbs distributed seeds.** For the target structure (((((...)))...(((...)))...), 1 000 sequences were generated, either uniformly or in the Boltzmann distribution, using IncaRNAtion [160].

* A diverse set of sequences allows for further filtering based on properties that are not easily necessarily imprinted at a computational level.

### 3.3.1 Sampling from the dual-partition function (117, 158, 160)

This last point is at the foundation of one of the PhD contribution of Vladimir Reinharz, supervised by Jérôme Waldispühl and myself [158, 160]. Namely, as illustrated in Figure 3.3 for a toy structure, the Boltzmann probability typically correlate strongly negatively meaning with the free-energy. This means that sequences with low free-energy (positive design) is also likely to optimize classic objectives of design, and thus represents a potentially promising approach for the inverse folding problem. This property was exploited by [28], which started its local search from the MFE sequence for a given structure, leading to improved performances.

We considered the dual-partition function $\mathcal{Z}_{S^*}$ of a structure $S^*$, defined as:

$$ \mathcal{Z}_{S^*} = \sum_{w \text{ compatible with. } S^*} e^{-\beta \cdot E(w, S^*)} $$

where $E(w, S)$ represents a simplified free-energy model, restricted to stacked base pairs, and $\beta$ is an arbitrary constant. We then turned to a random generation of sequences in the Boltzmann-Gibbs distribution over sequences, where any sequence $w$, compatible
with $S^*$, is emitted with probability

$$ \mathbb{P}(w | S^*) = \frac{e^{-\beta \cdot E(w, S^*)}}{Z_{S^*}}. $$

We refer the reader to Chapter 2.3.4 for the algorithmic aspects of computing the dual partition function, and performing the stochastic backtrack required for the sampling.

**Remark 3.3.1:** Note that, while it constitutes an approximation of the established Turner energy model [185], our stacking-based model correlates very well with the former, and virtually does not impact the above probability. In particular, the lack of support for multiple loops provably does not impact the distribution. Indeed, the contribution $E_M$ of a multiple loop within the Turner model depends on its structural characteristics, but remains independent of its precise nucleotide content. Including such a contribution to our energy model would therefore increase the Boltzmann factor of any sequence by a factor $e^{-\beta \cdot E_M}$, homogeneous across sequences. Meanwhile, the partition function would be inflated by the exact same factor $e^{-\beta \cdot E_M}$, so that the two terms would cancel in the emission probability, leaving it unchanged.

This approach generates sequences which not only have, unsurprisingly, lower free-energies than those produced by uniform sampling, but also much larger Boltzmann probabilities, as illustrated by the left panel of Figure 3.3. However, further analyses (Figure 3.3, right panel) reveal that this increased suitability to the objectives of design comes at the cost of a **drift of the GC-content towards high values** (typically 75-80%). This drift is problematic since as wild-type sequences within living organisms often present medium or low GC-content, presumably to offer better transcription rates and/or structural plasticity [54].

As a countermeasure, we introduced a **rejection strategy**, an instance of multidimensional Boltzmann sampling (see Section 2.3.4 and 4.2.2), to generate sequences at a constrained GC-content, introducing a weight $\chi_{gc}$ to induce a distribution

$$ \mathbb{P}_{\chi_{gc}}(w | S^*) \propto \chi_{gc}^{|w|_G + |w|_C} \times e^{-\beta \cdot E(w, S^*)}. $$

As shown in Figure 3.4, the drift of Boltzmann-generated sequences towards high GC contents can be counterbalanced, while preserving much better free-energies than induced by the Boltzmann distribution. Moreover, this strategy induces Boltzmann probabilities that are much higher than induced by a uniform sampling, as well as a much better proportion of sequences adopting the target structure as its MFE. In other words, this **positive design strategy** constitutes a promising approach, possibly in combination with further **local optimizations**, for negative design.
<figure><img src="image_35.png" /><figcaption>Figure 3.4: Influence of targeted GC-content on characteristics of seeds. For the structure of Figure 3.3, IncaRNAtion [160] is used to produce Boltzmann-weighted sequences, constrained to feature ~50% GC content (blue), comparable to the one induced by a uniform random generation (green). GC-constrained sequence achieve much better free-energies and Boltzmann probabilities than uniformly-distributed sequences, while</figcaption></figure>

### 3.3.2 Avoiding and forcing motifs (210)

This weighted sampling strategy can also be extended to accommodate further constraints, such as collections of **forbidden ($F$) and mandatory ($M$) motifs**.

**Problem 12 (Motifs-constrained dual partition function):**

**Input:** Two sets $F, M \subset \Sigma^*$, respectively representing **forbidden** and **mandatory** sequence motifs; target structure $S \in S_n$

**Output:** The motifs-constrained dual partition function:

$$Z_{S^*, F, M} = \sum_{\substack{w \text{ compatible with. } S^* \\ \text{s. t. } f \notin w, \forall f \in F \\ \text{and } m \in w, \forall m \in M}} e^{-\beta \cdot E(w, S^*)}$$

This problem can be solved by constructing the deterministic finite-state automaton $A$ that recognizes words featuring all elements of $M$, and none of $F$. $A$ is built from specialized versions $A_{M'}$ of the Aho-Corasick automaton [2], which we call **components**, for all sets $F \cup M' \subseteq M$. Components are then collated, by redirecting any terminal state $m \in M'$ of $A_{M'}$ towards a suitable state in $A_{M'\setminus\{m\}}$. Thus, starting from $A_M$ and reading a word $w$, one ends up in the component $A'_{M'}$, where $M \setminus M'$ is exactly the set
<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /><figcaption>Figure 3.5: Construction of the automaton for sequence-constrained positive design. From lists of forbidden ($\mathcal{F}$) and mandatory ($\mathcal{M}$) words (A), the Aho-Corasick automaton (B) can be built in time linear on the accumulated length of words. The full automaton (C) ensures the presence of all words in $\mathcal{M}$, by constructing components, i.e. AC automata for $\mathcal{F}$ and subsets $\mathcal{M}' \subseteq \mathcal{M}$ of yet-to-be-generated words. Transitions between components are defined to account for new elements of $\mathcal{F}$, and the only accepting states are those of the component $\mathcal{M}' \neq \emptyset$. Finally, any terminal state in (B) associated with a word in $\mathcal{F}$ is rerouted to a ground/garbage state $\perp$.</figcaption></figure>

of mandatory words in $w$. Additionally, all final states associated with motifs in $\mathcal{F}$ are rerouted to a non-accepting ground state. The whole process is further explained and illustrated by Figure 3.5.

The number of states $|\mathcal{A}|$ of $\mathcal{A}$ grows like $\Theta(2^{|\mathcal{M}|} \sum_{m \in \mathcal{F} \cup \mathcal{M}} |m|)$. The exponential dependency on $|\mathcal{M}|$ seems, in some sense, unavoidable. Indeed, Räihä and Ukkonen [162] have shown the NP-hardness of deciding whether or not all motifs of a collection $\mathcal{M}$ ($\mathcal{F} = \emptyset$) can be represented within a string of length $n$.

Then, the progression within $\mathcal{A}$ can be simulated within the dynamic programming scheme through explicit convolutions of the state space. We give an example of such a construct in the example below.

### Example 11: Dual-partition function and constraints

In a classic **base-pairing energy model**, assigning a free-energy contribution $E_{a,b}$ to a base pair involving nucleotides a and b, the dual partition function for a structure $S^*$ is computed for any interval $[i, j]$ as:

$$ \forall i, j \in [1, n] : Z_{i,j} = \begin{cases} 1 & \text{if } i > j \\ \sum_{a \in \{A,U,C,G\}} Z_{i+1,j} & \text{if } i \text{ unpaired in } S^* \\ \sum_{(a,b) \in \mathcal{B}} e^{-\beta \cdot E_{a,b}} \times Z_{i+1,k+1} \times Z_{k+1,j} & \text{if } i \text{ paired to } k \text{ in } S^* \end{cases} $$

and $Z_{S^*} := Z_{1,n}$. Moreover, starting from $[1, n]$, only a linear number of subintervals $[i, j] \subset [1, n]$
is reached by the computation, so the recurrence can be computed in linear-time using dynamic programming.

Then we consider the adjunction of sequence constraints, forcing the presence/absence of motifs $M/F$ modeled by an automaton $A$ with state space $Q$, transition function $\delta$ and final states $F \subset Q$. It can be achieved by convolving the dynamic-programming schemes and the automaton $A$, essentially simulating its transitions. For all pair of states $q, q' \in Q$ and indices $i, j \in [1, n]$, one has :

$$Z_{i,j}^{q \to q'} = \begin{cases} 0 & \text{if } i > j \text{ and } q \neq q' \\ 1 & \text{if } i > j \text{ and } q = q' \\ \sum_{a \in \{A, U, C, G\}} Z_{i+1, j}^{\delta(q, a) \to q'} & \text{if } i \text{ is unpaired in } S^* \\ \sum_{(a, b) \in B} \sum_{q'' \in Q} e^{-\beta \cdot E_{a, b}} \times Z_{i+1, k+1}^{\delta(q, a) \to q''} \times Z_{k+1, j}^{\delta(q'', b) \to q'} & \text{if } i \text{ is paired to } k \text{ in } S^* \end{cases}$$

where $Z_{i,j}^{q \to q'}$ represents the dual partition function, restricted to words $w$ that are: i) compatible with the interval $[i, j]$ in $S^*$; and ii) starting from a state $q$, one ends in $q'$ upon reading $w$ in $A$.

### 3.3.3 Multiple structures and extended design principles (88, 90)

In the previous section, the computation of the dual partition function was greatly helped by the tree-like nature of the secondary structure, inducing recurrences at the origin of our dynamic-programming schemes. The presence of complex networks of constraints, for instance induced by **multiple secondary structures** for multistable RNAs [65] or to prevent the extension of helices [91], breaks the substructure property induced by a single secondary structure, and impedes the design of efficient algorithms.

In collaboration with Stefan Hammer, Wei Wang and Sebastian Will [88, 90], we considered a generalized version of the dual partition function, where multiple target structures $S^*, S^* \dots S_m^*$ are jointly considered, and the energy of a sequence is additively defined on a collection of functions $(F_1, F_2, \dots, F_k)$, called **contribution**. Each contribution $F$ **ob-serves** a set $\text{dep}(F) \subset [1, n]$ of positions in the designed sequence, and returns a numerical value accordingly. In other words, each contribution is a function $F: \sum^{\text{dep}(F)} \to \mathbb{R}$, whose values accumulate to define the **energy** of any given sequence.

**Problem 13 (General dual partition function):**

**Input:** Length $n$; Collection of **contributions** $F := [F_1, F_2, \dots, F_k]$; Collection of **target structures** $S^* := [S_1^*, S_2^*, \dots, S_m^*]$.

**Output:** The generalized dual partition function:

$$Z_{S^*, F} = \sum_{\substack{w \in \Sigma^m \\ w \text{ comp. with. } S^*, \forall S^* \in S^*}} e^{-\beta \cdot F(w)}$$

where $F(w) = \sum_{i=1}^k F_i([w_p | p \in \text{dep}(F_i)])$.
### 3.3.3.1 Complexity aspects

Complexity-wise, the problem is generally NP-hard, as it specializes into a variant of negative design studied by, and proven NP-hard, by Hellmuth et al. [91]. Namely, the authors introduced **extended shapes**, possibly described as secondary structures whose base-pairs are partitioned into **regular** base pairs, which should be populated by valid pairs of nucleotides, and **pseudo-edges** which should not be pairable. Given several such extended shapes, the **Min-Edge-Deletion in Shapes (MES) problem** asks for the minimal number of (regular or pseudo) edges to delete before a satisfying sequence assignment can be found, i.e. such that nucleotides on regular edges (i, j) (resp. pseudo edges (ĩ, j̃)) forming valid pairs in $\mathcal{B}$ (resp. invalid pairs in $U := \{AA, CC, GG, UU\}$). This problem essentially constitutes a specialization of the above problem, obtained by: i) setting $S^* := \emptyset$; ii) Choosing suitable contributions $F_{(i,j)}(a, b) = -1_{(a,b)\in\mathcal{B}}$ and $F_{(i,j)}(a, b) = -1_{(a,b)\in U}$. In this setting, $Z_{\emptyset,\mathcal{F}} = \sum_w e^{-\beta \cdot \#\text{DelEdges}(w)}$, and the existence of a sequence violating k edges can be decided from the dual partition function in polynomial time (e.g. through interpolation, or increasing $\beta$ so that co-optimal overwhelm the partition function).

Moreover, in a *uniform setting* ($\mathcal{F} := \emptyset$), a key object is the **dependency graph** $G_S$ of a collection of structures $S^* := [S_1^*, \dots, S_m^*]$ of length n, defined as

$$G_S := ([1, n], \cup_{i=1}^m S_i).$$

This graph is simply the union of the graphs of secondary structures, and must be **bipartite** for the collection to admit any design [73].

**Theorem 5 (Independent set and multidesign):** Given several secondary structures $S^* := [S_1^*, \dots, S_m^*]$ such that $G_S$ is **bipartite**, the **dual partition function** $Z_{S^*,\emptyset}$ is exactly twice the number of independent sets in $G_S$.

It follows that computing the general partition function is as hard as **counting the number of independent sets in a bipartite graph (#BIS)**, a thoroughly-studied #P-hard problem [79].

Finally, we described in Hammer et al. [88, 90] a general **Fixed Parameter Tractable** solution. The entire network of contributions is jointly considered, and automatically transformed into a set of DP equations through a **tree decomposition**. Resulting DP schemes contrast with most of the schemes described in this document due to the automated nature of their design, but also due to their optimality in some well-defined sense.
# Chapter 4

## Constrained random generation through rejection and beyond

Generic approaches for the random generation of combinatorial objects rely on the existence of a specification for the combinatorial class. However, the design of a specification for a given combinatorial class can be a highly complex task, or altogether be provably impossible for a given a restricted set of constructors. In the case of combinatorial classes that are naturally expressible as languages, one can use classic pumping-lemma to rule out the existence of regular or context-free specifications. However, such results do not preclude the existence of a grammar for a language, which itself would be in bijection with the language of interest, as shown in the following example.

**Example 12: Copy language**

Consider the 3-copy language

$$L_{3C} = \{a^n.b^n.c^n \mid n \in \mathbb{N}\}$$

A textbook application of the pumping lemma (algebraic version) shows that the language is not context-free: there is no context-free grammar that generates $L_{3C}$. However, the language

$$L_3 = \{(a.b.c)^n \mid n \in \mathbb{N}\}$$

can be generated using a context-free grammar (or even a simpler finite-state automaton), and admits a trivial bijection with $L_{3C}$. Such a bijection may simply push occurrences of a (resp. b and c) to the beginning (resp. ending and middle), preserving the relative order within occurrences of a given letter. This transform can is easily reversible, and can be implemented in linear time.

Combinatorial properties of languages can be used to rule out any such workaround and, more generally, establish the unsuitability of classes of grammars for the generation of some languages.

For instance, Flajolet [67] showed that some context-free language are **intrinsically ambiguous**, i.e. any grammar that generates such a language must be ambiguous. Indeed, languages generated from a context-free grammar are associated with ordinary generating functions that are algebraic. However, certain context-free languages are known
to admit transcendental generating functions, and thus cannot admit an unambiguous context-free grammar, as shown in Example 13.

**Example 13: A simple intrinsically ambiguous language<sup>a</sup>**

Consider the language

$$L_{IA} = \{w \in \{a, b, c\}^* \mid |w|_a \neq |w|_b \text{ or } |w|_b \neq |w|_c\}.$$

It is generated by a context-free grammar (axiom S) $S \to T_{a<b} | T_{b<a} | T_{b<c} | T_{c<b}$, which breaks down the language into four (overlapping, thus inducing an ambiguous grammar) sub languages that are context free, and whose union covers all possible cases. For instance, the non-terminal $T_{a<b}$ is tasked with generating the language

$$L_{a<b} = \{w \in \{a, b, c\}^* \mid |w|_a < |w|_b\}$$

using the following set of production rules

$$T_{a<b} \to a T_{a=b} b T_{a<b} \mid b T_{a=b} a T_{a<b} \mid b T_{a<b} \mid b T_{a=b} \mid c T_{a<b}$$

$$T_{a=b} \to a T_{a=b} b T_{a=b} \mid b T_{a=b} a T_{a=b} \mid c T_{a=b} \mid \varepsilon$$

which essentially amounts to a case distinction based on the first encountered character, reminiscent of a *first passage decomposition*. Similar productions rules can be obtained for the languages generated by $T_{b<a}$, $T_{b<c}$ and $T_{c<b}$ respectively, establishing the context-free nature of $L_{IA}$.

To study the generating function of $L_{IA}$, remark that $L_{IA} = \{a, b, c\}^* \setminus L_{\overline{IA}}$, where

$$L_{\overline{IA}} = \{w \in \{a, b, c\}^* \mid |w|_a = |w|_b = |w|_c\}$$

The ordinary generating functions $S_{IA}(z)$ and $S_{\overline{IA}}(z)$ are thus related through

$$S_{IA}(z) = \frac{1}{1 - 3z} - S_{\overline{IA}}(z).$$

Remind that algebraic generating functions are closed under addition/subtraction, so $S_{IA}(z)$ is algebraic only if $S_{\overline{IA}}(z)$ is algebraic itself. However, $S_{\overline{IA}}(z)$ can be expressed explicitly as

$$S_{\overline{IA}}(z) = \sum_{n \ge 0} \binom{3n}{n, n, n} z^{3n}.$$

It follows from Stirling's expansion that, for $n \equiv 0 \mod 3$, one has

$$[z^n] S_{\overline{IA}}(z) = \frac{n!}{(n/3)! \times (n/3)! \times (n/3)!} = \frac{3 \times 3^n}{2\pi \times n} (1 + O(1/n)).$$

Such an asymptotic expansion is not compatible with the universal asymptotic form of coefficients of rational/algebraic generating functions, as seen in the Flajolet/Sedgewick *bible* [69], or a recently refined result by Banderier and Drmota [11] (see Remark 5.1.2).

$S_{\overline{IA}}(z)$ is therefore not algebraic, and neither is $S_{IA}(z)$. Since generating functions associated with unambiguous context-free grammars are either rational or algebraic, we conclude that no such grammar exists, and $L_{IA}$ is context-free, yet *intrinsically ambiguous*.

More generally, combinatorial classes associated with transcendental generating functions, regardless of their language-theoretic properties, cannot be modeled using the standard operators introduced by the symbolic method. In this chapter, we remind,

<footnote><sup>a</sup>Example borrowed from the *magisterial exposition* of Cyril Nicaud.</footnote>
generalize and exploit versatile rejection principles for the generation of constrained
combinatorial objects which are hard to specify (at least in their constrained version).

**Outline.** Section 4.1 starts by reminding the general objectives and principles of rejection methods for random generation. Section 4.2 describes the main principles behind the general technique called multidimensional Boltzmann sampling, whose application to selected problems in Bioinformatics is described in Section 4.3. Section 4.4 provides an analysis of the level of redundancy within weighted samples, and describes possible countermeasures to such an unproductive phenomenon.

The following summarizes, and attempts to unify, contributions described within the
following list of articles published in journals and/or presented at conferences.

Associated contributions

M. Bousquet-Mélou and Y. Ponty. Culminating paths. Discrete Mathematics and Theoretical Computer Science, 10(2):125-152, 2008

A. Denise, Y. Ponty, and M. Termier. Controlled non uniform random generation of decomposable structures.
**Theoretical Computer Science**, 411(40-42):3527-3552, 2010

O. Bodini and Y. Ponty. Multi-dimensional Boltzmann Sampling of Languages. In **AOFA 2010**, volume AM of *DMTCS Proceedings*, pages 49–64, Vienna, Austria, June 2010. Discrete Mathematics and Theoretical Computer Science

D. Gardy and Y. Ponty. Weighted random generation of context-free languages: Analysis of collisions in random urn occupancy models. In **GASCOM 2010**, page 14pp, Montréal, Canada, Sept. 2010. LACIM, UQAM

J. Waldispühl and Y. Ponty. An unbiased adaptive sampling algorithm for the exploration of RNA mutational landscapes under evolutionary pressure. Journal of Computational Biology, 18(11):1465–79, Nov. 2011

<figure><img src="image_39.png" /></figure>

J. Waldispühl and Y. Ponty. An unbiased adaptive sampling algorithm for the exploration of RNA mutational landscapes under evolutionary pressure. In *RECOMB 2011*, volume 6577 of *Lecture Notes in Computer Science*, pages 501–515, Vancouver, Canada, Mar. 2011. Springer Berlin / Heidelberg

J. Du Boisberranger, D. Gardy, and Y. Ponty. The weighted words collector. In **AOFA 2012**, volume AQ, pages 243–264, Montreal, Canada, June 2012. DMTCS

C. Banderier, O. Bodini, Y. Ponty, and H. Tafat. On the diversity of pattern distributions in rational language. In **ANALCO 2012**, pages 107–116, Kyoto, Japan, Jan. 2012. Omnipress

A. Lorenz and Y. Ponty. Non-redundant random generation algorithms for weighted context-free languages.
*Theoretical Computer Science*, 502:177–194, Sept. 2013

<figure><img src="image_40.png" /></figure>

Y. Ponty. Non-redundant random generation from weighted context-free languages. In **GASCOM 2008**, page 12pp, Bibbiena, Italy, 2008

V. Reinharz, Y. Ponty, and J. Waldispühl. A weighted sampling algorithm for the design of RNA sequences with targeted secondary structure and nucleotide distribution. *Bioinformatics*, 29(13):i308–15, July 2013

<figure><img src="image_41.png" /></figure>

V. Reinharz, Y. Ponty, and J. Waldispühl. A weighted sampling algorithm for the design of RNA sequences with targeted secondary structure and nucleotides distribution. In **ISMB/ECCB 2013**, Berlin, Germany, July 2013

Y. Zhang, Y. Ponty, M. Blanchette, E. Lecuyer, and J. Waldispühl. SPARCS: a web server to analyze (un)structured regions in coding RNA sequences. Nucleic Acids Research, 41(Web Server issue):W480–5, July 2013

J. Lumbroso, M. Mishna, and Y. Ponty. Taming reluctant random walks in the positive quadrant. In **GASCOM 2016**, volume 59 of *Electronic Notes in Discrete Mathematics*, pages 99 – 114, Bastia, France, 2017
J. Michálik, H. Touzet, and Y. Ponty. Efficient approximations of RNA kinetics landscape using non-redundant sampling. In ISMB/ECCB 2017, volume 33, pages i283 – i292, Prague, Czech Republic, 2017

S. Hammer, W. Wang, S. Will, and Y. Ponty. Fixed-parameter tractable sampling for RNA design with multiple target structures. BMC Bioinformatics, 20(1):209, Dec. 2019

S. Hammer, Y. Ponty, W. Wang, and S. Will. Fixed-Parameter Tractable Sampling for RNA Design with Multiple Target Structures. In RECOMB 2018, Paris, France, 2018

## 4.1 It's not you, it's me: The unreasonable power of rejection

Before we move any further, let us define precisely the problem of random generation. We are tasked with the random generation of objects of size n within a challenging class A (i.e. hard to specify, e.g. using a grammar), according to a predefined probability distribution $p_w^*$.

**Problem 14 (Random generation for A):**

**Input:** Number m of objects

**Output:** Sequence of $[w_1, \dots, w_m] \in A^m$ such that

$$ \mathbb{P}_A(W = w) = p_w^* $$

The provable inexistence of a specification can be frequently circumvented at a reasonable cost through a **rejection strategy**. Its main idea is to consider another simpler set B, chosen to enforce the following properties:

1. **Correctness:** $A \subseteq B$ (or $A \subseteq B'$, admitting a computable bijection with B);

2. **Effectiveness:** Objects from B can be generated (uniformly) at random;

3. **Efficiency:** The distribution over B is reasonably close to its target over A.

In this setting, the rejection algorithm consists in an iterated generation of objects from a superset $B \supseteq A$, discarding any object in $B \setminus A$, until an object of A is obtained and returned.

The subsequent rejections act as a conditioning, and the **emission probability** $p_w$ of a given object $w \in A$ simply becomes

$$ p_w = \mathbb{P}_B(W = w | W \in A) = \frac{\mathbb{P}_B(W = w)}{\mathbb{P}_B(W \in A)} = \frac{\mathbb{P}_B(W = w)}{\sum_{w' \in A} \mathbb{P}_B(W = w')} $$

where $\mathbb{P}_B(W = w)$ indicates the probability of generating w using the above-mentioned generator for B.
Typically, when a **uniform distribution** $p_w^* = 1/|A|$ is targeted, and enforced over $B$ by the initial generator, one has $P_B(W = w) = 1/|B|$, and the emission probability simplifies into

$$p_w = \frac{1/|B|}{|A|/|B|} = \frac{1}{|A|} = p_w^*.$$

In other words, a uniform generator for $B$ readily yields a uniform generator for $A$.

This correctness generalizes **beyond uniform distributions**, to any case where the generated and targeted probability distributions are proportional over $A$:

$$\exists \kappa \in \mathbb{R} \text{ such that } \forall w \in A, \quad P_B(W = w) = \frac{p_w^*}{\kappa} \qquad (4.1)$$

where $p_w^*$ denotes the targeted probability distribution. Again, the emission probability simplifies into

$$p_w = \frac{p_w^*/\kappa}{\sum_{w' \in A} p_w^*/\kappa} = p_w^*,$$

and the rejection generator enforces the targeted probability distribution.

### 4.1.1 Correcting a (bounded) bias

A similar strategy can be used when the distribution of $B$, renormalized over $A$ (i.e. after rejection), does not fully coincide with the targeted distribution:

$$\exists w \in A \text{ such that } p_w \neq p_w^*.$$

Denote as $w_\Theta$ the object in $A$ which is most unfavorably discriminated against by the initial generation, such that

$$w_\Theta := \underset{w \in A}{\operatorname{argmin}} \delta(w) \quad \text{with} \quad \delta(w) := \frac{p_w}{p_w^*}.$$

Then, following the ideas of [142], the rejection algorithm is coupled with an additional rejection step, accepting the generated $w \in A$ with probability

$$q_w = \frac{\delta(w_\Theta)}{\delta(w)} \le 1 \qquad (4.2)$$

and rejecting it otherwise. The emission probability $p'_w$ of the combined rejection is then

$$p'_w = \frac{p_w \times q_w}{\sum_{w' \in A} p_{w'} \times q_{w'}} = \frac{p_w^*}{\sum_{w' \in A} p_{w'}^*} = p_w^*$$

so that any generator for $B$ such that $P_B(W = w) > 0, \forall w \in A$, can be used (potentially at the cost of frequent rejections) to generate from $A$ within the targeted distribution.
### 4.1.2 Complexity aspects

Given the existence of unbounded sequences of rejections the worst-case complexity of rejection generators is infinite. However, such scenarios are associated with abysmally small probabilities and, overall, populate a negligible proportion of the probability space, so worst-case analyses provides little insight into the concrete requirements of rejection-based methods. It is thus arguably more relevant to consider the trade-off between cost and probability, and **average-case complexities** are thus usually sought to measure the practical complexity of generators.

The main contributing factor to the average-case complexity usually resides in the repeated (discarded) calls to the generator for $\mathcal{B}$. Denote by $T$ the **number of attempts**, i.e. the number of time the generator for $\mathcal{B}$ must be called before generating an object from the targeted distribution over $\mathcal{A}$, then we **generally** have

$$ \mathbb{E}(T) = \frac{1}{\sum_{w'' \in \mathcal{A}} \mathbb{P}_{\mathcal{B}}(W = w'') \times q_{w''}} \\ \le \frac{1}{\mathbb{P}_{\mathcal{B}}(W \in \mathcal{A})} \times \frac{p_{w_{\ominus}}^{*}}{\mathbb{P}_{\mathcal{B}}(W = w_{\ominus})} \times \frac{p_{w_{\oplus}}^{*}}{\mathbb{P}_{\mathcal{B}}(W = w_{\oplus})}, $$

with $w_{\oplus} = \max_{w \in \mathcal{A}} \delta(w)$ being the object whose emission probability is the most positively biased in comparison to the targeted distribution.

Note that, when the **uniform distribution** is targeted over $\mathcal{A}$, and exactly enforced over $\mathcal{B}$, we have $q_w = 1, \forall w \in \mathcal{A}$, and the expected number of attempts greatly simplifies into

$$ \mathbb{E}(T) = \frac{|\mathcal{B}|}{|\mathcal{A}|}. $$

#### Example 14: 2D walks with small steps (122)

Many types of 2D walks are associated with non-algebraic generating functions when confined to the positive quadrant.

For instance, Mishna and Rechnitzer [137] showed that the generating function of positive 2D walks taking steps in $S = \{\llcorner, \rfloor, \rfloor\}$ is not D-finite and, in particular, non-algebraic. However, the number $w_n$ of such walks taking n steps admits an asymptotic equivalent in

$$ c_n = \kappa 3^n + \mathcal{O}(\sqrt{8^n}) $$

where $\kappa = 0.173...$ is a computable constant.

Now consider the set of unconstrained walks over $\{\llcorner, \rfloor, \rfloor\}^n$, consisting of exactly $3^n$ objects. A trivial algorithm for the uniform generation simply chooses n steps, each one uniformly at random.

Clearly, the set of unconstrained walks (set $\mathcal{B}$) includes the set of positive 2D walks (set $\mathcal{A}$). A rejection algorithm can then be envisioned, generating unconstrained walks from $\mathcal{B}$ until a positive walks from $\mathcal{A}$ is found. The expected number of attempts from this rejection algorithm is then asymptotically