equivalent to

$$E(T) = \frac{|\mathcal{B}|}{|\mathcal{A}|} = \frac{3^n}{\kappa 3^n} \approx 5.77...$$

an, in particular, does not grow with n, leading to a linear-time rejection algorithm.

With Marni Mishna and Jérémie Lumbroso, we generalized this strategy [122] to more complex small steps set $S \subset [-1, 0, 1] \times [-1, 0, 1]$, crucially allowing for an efficient generation in cases of **reluctant walks**, i.e. walks having a natural attraction, called **drift**, towards the axes. This tendency to violate the positivity condition in turn leads to a number of positive walks exponentially smaller than $|S|^n$, thus rejecting from unconstrained walks would yield an exponential expected time algorithm.

However, it as been shown in Duraj [61] that the number $w_n$ of positive 2D walks is such that $w_n \in \Theta(\alpha^n n^\beta)$, for two constants $\alpha < |S|$ and $\beta$ that only depend on the step set $S$. Moreover, this exponential growth constant $\alpha$ is shared by a well-chosen half-plane, including the positive quarter-plane, whose slope m can be determined numerically (or even exactly). When such a slope m is rational, then the projection of the steps onto the chosen-plane is a rational number, and efficient random generation algorithms for 1D positive walks can be used to generate walks in the half plane. Those walks are then lifted onto the 2D plane, and an additional rejection step ensures positivity at the 2D level.

Since the exponential growth of the 1D and 2D walks coincide, then the rejection only contributes an additional polynomial factor. In other words, the algorithm that generates 1D positive walks and rejects those that escape the 2D positive quadrant, has expected time per generation bounded by a polynomial in n.

Finally, we consider the **computational complexity**, expressed in an arbitrary unit (#arithmetic operations, bit-complexity [49]), and denoted as a random variable C.

In simple – yet frequent – cases, objects of a fixed size n are generated, and the time complexity $c(w)$ required by the generation of an object can be seen as deterministic and only depending on n. Denoting by $c(n)$ the complexity of generating objects in $\mathcal{B}$, one gets

$$E(C) = \frac{c(n)}{\mathbb{P}_{\mathcal{B}}(W \in \mathcal{A})} \qquad (4.3)$$

More generally, the average-case complexity of the rejection algorithm follows

$$E(C) = \frac{\mathbb{E}_{\mathcal{B}}(c(W))}{\mathbb{P}_{\mathcal{B}}(W \in \mathcal{A})} \qquad (4.4)$$

where W is a random object in $\mathcal{B}$. This formula also holds when the expected cost of generating objects in $\mathcal{A}$ and $\mathcal{B} \setminus \mathcal{A}$ differ, owing to the additivity, in expectation, of (dependent) random variables.

### Example 15: Anticipated rejection/Florentine algorithm for Motzkin walks (13)

An example of such an analysis can be found in the generation of **prefixes of Motzkin walks**, i.e. random positive paths using ↗, → and ↘ steps, associate to 2D vectors (+1, +1), (+1, 0) and (+1, -1) respectively. The textbook algorithm, proposed by Barcucci, Pinzani and Sprugnoli [13], uses an efficient **anticipated rejection** scheme. Its algorithmic part is straightforward, and simply consists, starting from 0, of choosing ↑, → or ↓ with equal probability at each stage, rejecting the prefix as soon
as a negative ordinates is reached, until a requested length n is reached.

The analysis of its average-case complexity is slightly more involved. The expected number of attempts depends on the asymptotics of the numbers of positive paths ($A_n$) and unconstrained paths ($B_n$) respectively. While we clearly have $|B_n| = 3^n$, it is a folklore result that $|A_n| \in \Theta(\frac{3^n}{\sqrt{n}})$, so the probability of accepting an object is

$$P_B(W \in A) = \frac{|A|_n}{|B|_n} \in \Theta(1/\sqrt{n}).$$

The expected number of attempts thus grows in $\Theta(\sqrt{n})$ and a classic rejection scheme, generating a n step unconstrained walk followed by a positivity test, would yield and overall complexity in $\Theta(n\sqrt{n})$.

Anticipated rejection has a sizable impact on the expected generation cost, but its analysis cost is slightly more involved. Denote by $b_n$ the number of positive walks of length n, and by $d_n$ the number of positive walks ending in (n, 0). Being rejected before reaching n steps means having a positive prefix ending at (k, 0) of size k < n steps, followed by a ↘ step, with a cost of k rejections being incurred by such a rejection. The expected cost of a generation is therefore given by

$$E(c(W)) = \sum_{k<n} k \frac{d_{k-1}}{3k} + n \frac{b_n}{3n},$$

where the first term denotes the expected cost of (anticipated) rejections, and the second term denoting the contribution of positive paths. It then follows from

$$d_n \in \Theta\left(\frac{3^n}{n\sqrt{n}}\right) \quad \text{and} \quad b_n \in \Theta\left(\frac{3^n}{\sqrt{n}}\right),$$

that the first term is bounded by a constant, while the second grows like $\Theta(\sqrt{n})$, and thus

$$E(c(W)) \in \Theta(\sqrt{n}).$$

A direct application of (4.4) allows us to conclude that the overall average-case complexity is $\Theta(n)$. This result is optimal (up to constants) since the number of bits required despite an expected $\Theta(\sqrt{n})$ number of attempts.

A similar reasoning, and complexity speedup, more generally holds for a very large class of combinatorial structures obtained by step-wise extensions, such that the extension of a rejected object is itself rejected, as shown by Denise [48].

### Example 16: Culminating paths (21)

Using anticipated rejection, to spare the cost of rejecting a fully-generated object, was crucial in the design, with Mireille Bousquet-Mélou, of an efficient algorithm for the uniform generation of **culminating paths** [21]. A culminating path is a unidimensional walk which takes steps +1 and -1, starts at position 0, and must both: 1) remain positive; and 2) be culminating, i.e. reaches its maximum height on its final step.

Now consider the **mirror image** of a positive walk, obtained by reversing its sequence of steps, and note that it is now culminating (but may not be positive). It follows that any culminating path $w$ of length $n$ can be factored into $w = w^{\perp} \cdot w^{\top}$, where $w^{\perp}$ is a positive path of size $\lceil n/2 \rceil$, and $w^{\top}$ is the mirror image of some positive path of length $\lfloor n/2 \rfloor$. In other words, the set $B_n := W_{\lceil n/2 \rceil}^{\perp} \times W_{\lfloor n/2 \rfloor}^{\top}$ where $W_n^{\perp}$ (resp. $W_n^{\top}$) is the set of positive (resp. culminating) walks of length $n$, represents a superset of $A_n$ the set of culminating walks of length $n$.
Using the kernel method, we could establish that $|A_n| \in \Theta(2^n/n)$ [21], and it is a folklore result that $|W_m^\perp| = |W_m^T| \in \Theta(2^n/\sqrt{n})$, so that the expected number of attempts is asymptotically

$$E(T) = \frac{|W_{[n/2]}^\perp| \times |W_{[n/2]}^T|}{|A_n|} \in \Theta(1).$$

Moreover, the generation of a positive walk, or its reversal can be done in linear expected time using anticipated rejection. Overall, we get that the simple algorithm that generates two half walks, reverses the second one, and glues them together, restarting until the result is a culminating path, has linear average-case complexity.

## 4.2 Rejecting in dimension one and then some...

### 4.2.1 (Combinatorial) Boltzmann sampling

**Boltzmann sampling** [60] can be seen as a special instance of rejection sampling. It revisits the random generation of combinatorial objects of a given size n described by a specification B. Instead of precomputing the number of objects accessible through local alternatives using costly convolution products, it initially relaxes the size constraint, drawing objects of any size in B until an object of suitable size is found.

More precisely, a **Boltzmann generator** emits an object $\omega \in B$ with probability

$$P_x(W = \omega) = \frac{\chi^{|\omega|}}{B(\chi)} \quad (4.5)$$

where $B(z)$ is the generating function of B and $\chi$ is a real-valued **Boltzmann parameter**. Such a distribution can be seen as a special instance of the Boltzmann distribution in statistical mechanics, where the free-energy of an object $\omega \in B$ is arbitrarily set to $E(\omega) := -kT \log x|\omega|$. This leads to associating a Boltzmann factor of $e^{E(\omega)/kT} = \chi^{|\omega|}$ to each object, in which one recognizes the distribution of Equation (4.5).

A careful tuning of $\chi$ ensures a *reasonable* probability of acceptance and, in turn, an average-case complexity scaling in $O(n^2)$ for generating an object of size n (generating an object of approximate size in $[(1 - \varepsilon) \cdot n, (1 + \varepsilon) \cdot n]$ can be performed in $\Theta(n)$ time).

The appeal of Boltzmann samplers to a large subset of the Analysis of Algorithms (AOFA) community, can be attributed to two main factors:

* The elegance, and relative simplicity, of its algorithmic aspects: As shown in Table 4.1, generators for Cartesian products are deceptively simple to express. Their implementation, however, requires an evaluation of the **Boltzmann Oracle**, the value of the generating functions. While this task is relatively trivial when a closed form (or, generally, computable) formula is known for the generating function, its
<table><caption>Table 4.1: Uniform Boltzmann generators for context-free languages.</caption><thead><tr><th>Rule type</th><th>Generating function</th><th>Boltzmann generator ΓV<sub>x</sub></th></tr></thead><tbody><tr><td>Product v → a.b</td><td>V(z) = A(z) × B(z)</td><td><b>Return</b> ΓA<sub>x</sub>.ΓB<sub>x</sub></td></tr><tr><td>Union v → a, v → b</td><td>V(z) = A(z) + B(z)</td><td><b>If</b> Bernoulli(A(x)/V(x)) : <b>Return</b> ΓA<sub>x</sub><br><b>Else</b> : <b>Return</b> ΓB<sub>x</sub></td></tr><tr><td>Terminal v → t</td><td>V(z) = z</td><td><b>Return</b> t</td></tr><tr><td>Epsilon v → ε</td><td>V(z) = 1</td><td><b>Return</b> ε</td></tr></tbody></table>

automation for general specifications is more involved, yet possible using a mixture of formal calculus and numerical iterations [145]). The optimization of the Boltzmann parameter can also be addressed in its full generality using efficient numerical iterations[144];

* Its strong connection with analytic combinatorics: The time complexity of Boltzmann samplers is strongly related to the distribution of lengths under the Boltzmann distribution. The design of Boltzmann samplers thus motivates extensions of the symbolic method to complex operators (e.g. box-product, shuffle product...) that lead to generic random generation principles, and to deeper characterizations of the rich singularities that arise from those operators.

## 4.2.2 Multidimensional Boltzmann sampling

In collaboration with Olivier Bodini [18], we introduced<fnref n="1" /> and studied an extension of Boltzmann sampling, dubbed **multi-dimensional Boltzmann sampling**. This technique addresses cases where one wishes to express cardinality constraints on distinguished atoms/letters (or combination thereof). For instance, given a language L over an alphabet $Σ := \{a, b, c\}$, one may be interested in the subset of L where the number occurrences of a, b and c are equal, reducible to a pair of constraints $\{a = n/3; b = n/3\}$. Enforcing such constraints generally escapes the expressiveness of classic operators of the symbolic method, preventing the definition of a specification and, in turn, the use of generic generators.

### 4.2.2.1 General principle

The idea of generalizing the Boltzmann rejection scheme to higher dimensions stems from the observation that additive parameters (e.g. used to capture the #occurrence of atoms/letters, #use of production rules...) occurring in combinatorial specifications tend to asymptotically follow normal distributions. Moreover, those distributions have

<footnote n="1">¹This natural extension was already discussed in the conclusion of the seminal Boltzmann paper [60], yet remained unexplored for almost a decade prior to our work.</footnote>
linear mean and variance. Such distributions are typically robust to the introduction of **weights**, also called **parameters** in subsequent work, which can be used to **shift the means**, thus providing a probabilistic control over the average composition of generated objects.

Once fitted to match desired expectations, a rejection algorithm performs a weighted random generation, only accepting objects with suitable composition. The concentration of the distribution usually induces a polynomial expected number of attempts, typically growing in $Θ(n^{k/2})$ when k is the number of cardinality constraints. Let us illustrate our general approach with a brief introductory example.

### Example 17: Multidimensional Boltzmann sampler for balanced Motzkin language

Consider the Motzkin language $M \subseteq \{(\cdot, \cdot), \bullet\}^*$ of well-parenthesized expressions interleaved with dots, formally defined as:

$$M = \{\omega \in \{(\cdot, \cdot), \bullet\}^* \mid |\omega|_{\cdot} = |\omega|_{\bullet} \text{ and } |\omega'|_{\cdot} \ge |\omega'|_{\bullet}, \forall \omega' \sqsubseteq \omega \text{ (prefix of } \omega)\}.$$

This language is context-free, and generated by the grammar:

$$S \to (S) S \mid \bullet S \mid \varepsilon$$

Now, consider the balanced Motzkin language $M_{2B}$ the subset of the Motzkin words featuring exactly 50% of dots (•), i.e.

$$M_{2B} = \{ \omega \in M \mid (|\omega| \equiv 0 \mod 2) \text{ and } |\omega|_{\bullet} = \frac{|\omega|}{2} \}$$

Clearly, $M_{2B}$ is not context-free. Moreover, the targeted 1/2 proportion of dots (•) characters differs substantially from the 1/3 proportion expected in large uniformly-distributed Motzkin words [50]. It follows from Drmota's theorem [56] that a naive rejection algorithm, generating words of length n from M until a word in $M_{2B}$ is obtained, would require a number of attempts growing exponentially with n.

However, as previously described [50], a weight $\pi \in \mathbb{R}^+$ can be associated to the dot (•), to induce a **weighted distribution** onto M, in which the emission probabilities become

$$P_{\pi}(W = \omega) = \frac{\pi^{|\omega|\cdot}}{\sum_{\omega' \in M} \pi^{|\omega'|\cdot}} \quad (4.6)$$

Generating within such a distribution can be performed in time essentially $O(n \log n)$, after a precomputation in $O(n^{1+o(1)})$, using generic generators [50].

The expected proportion of dots then becomes a function of $\pi$, allowing the user to **shift the mean** by changing the value of $\pi$. Namely, the asymptotic proportion of dots is given by

$$E_{\pi}(|W|_{\bullet}) = \frac{\pi}{2 + \pi} (1 + O(1/\sqrt{n}))$$

and thus achieves the targeted 1/2 for $\pi = 2$. More importantly, if follows from the multidimensional version of Drmota's theorem [56] that, within the weighted version of M, the number of dots asymptotically follows a normal distribution with mean n/2, variance $\sigma n$, for some positive constant $\sigma$. The probability of generating a word in $M_{2B}$ is thus in $Θ(1/\sqrt{n})$.

Now, consider the rejection algorithm which repeatedly generates words in M of length n from the weighted distribution (4.6), setting $\pi = 2$, until a word in $M_{2B}$ is obtained. This algorithm is uniform
within the restriction of $M_{2B}$ to words of size n. Moreover, its expected number of attempts is in $\Theta(\sqrt{n})$, leading to an overall complexity in $\Theta(n\sqrt{n} \log n)$.

### 4.2.2.2 Problem statement and algorithm

The problem addressed by multidimensional Boltzmann sampling is the uniform random generation of words/objects of size n satisfying:

* Specification/grammar C for a combinatorial class over atoms $t_1, t_2...$ (a.k.a. letters, or terminal symbols);

* Composition constraints of the form

$$ \left\{ \frac{|W|_{t_i}}{n} = r_i \in [0, 1] \right\}_{i=1}^{k} \qquad (4.7) $$

for a subset of k atoms, denoted {$t_1, ..., t_k$} for the sake of convenience.

A multidimensional Boltzmann sampler [18] then proceeds in three stages:

1. Calibration of weights/parameters $\pi_1, \pi_2...$, chosen such that

$$ \left\{ \mu_i := \frac{\mathbb{E}_{\pi_1, \pi_2, \dots} (|W|_{t_i})}{n} = r_i \in [0, 1] \right\}_{i=1}^{k} \qquad (4.8) $$

2. Random generation of objects from C within the $\pi$-weighted distribution

$$ \mathbb{P}_{\pi_1, \pi_2, \dots} (W = \omega) = \frac{\prod_{i=1}^{k} \pi_i^{\omega|z_i}}{[z^n] C(z; \pi_1, \pi_2, \dots)} \qquad (4.9) $$

where $C(z; u_1, u_2...)$ represents the **multivariate (weighted) generating function** marking the occurrences of the atoms

$$ C_{\pi}(z; u_1, u_2, \dots) = \sum_{w \in C} z^{|w|} \prod_{i=1}^{k} \pi_i^{|w|_{t_i}} u_i^{|w|_{t_i}} \qquad (4.10) $$

3. Rejection of objects having unsuitable composition (and size) until adequate ones are produced.

The requirement of an exact size n for the generated objects can be relaxed to accept words having more diverse lengths in $[(1 - \epsilon) \cdot n, (1 + \epsilon) \cdot n]$ for some constant $\epsilon > 0$. Similarly, words whose composition departs slightly from the expressed constraints, such that $\#t_i/n \in [r_i - \epsilon', r_i + \epsilon']$, can be deemed acceptable.
**Remark 4.2.1:** Note that these relaxed statements only enforce uniformity within a given size/composition class, and do not constrain the respective probabilities of two objects having different length/composition. Such a slight bias is generally deemed acceptable in contexts where statistical estimates are computed from generated structures.

Moreover those emission probabilities could, in principle, be corrected by a subsequent rejection to ensure an overall uniform distribution, as described above in Equation (4.2). However, this additional rejection would increase the expected complexity by an exponential factor, roughly equal to $(\prod_i \max(\pi_i, 1/\pi_i))^{2\epsilon n}$.

**Weight calibration.** This step can be performed using two distinct approaches:

* Symbolically, it is often possible to find a simple expression relating the asymptotic proportion of each atom to the weight. For instance, in context-free languages generated by strongly-connected grammars, the dominant singularity is always of *square-root* type, and the dominant singularity $\rho_{\pi_1, \pi_2}...$ can be expressed as a function of the weights. The expected composition is then given by

$$ \lim_{n \to \infty} \mathbb{E} (|W|_{t_i}/n \mid |W| = n) = -\pi_i \frac{\frac{\partial}{\partial \pi_i} \rho_{\pi_1, \pi_2} \dots}{\rho_{\pi_1, \pi_2} \dots}, \forall i \in [1, k]. $$

This system can then be solved to obtain a closed-form formula, expressing the solution weights as a function of the targeted proportions, as shown in Example 17 and [50].

* Numerically, an iteration was initially proposed [18], whose efficient convergence was conditioned by the choice of a suitable initial weight vector.

This work was later revisited by Bendkowski, Bodini and Dovgal [15], leading to a more general and efficient solution. In their work, the authors were able to reformulate the calibration as a convex optimization problem, solved using the interior-point method of Nesterov and Nemirovskii [140] in time polynomial to the inverse of the required precision. The authors completed their study with estimates for the precision required to ensure a good acceptance probability, leading to an overall polynomial tuning of multidimensional Boltzmann sampling.

Whenever feasible, a symbolic approach is arguably more elegant and efficient, leading to a trivial algorithm for computing the weights. Unfortunately, it is easy to find examples where such a strategy involves solving systems whose complexity lies beyond the current capacities of symbolic mathematical engines (e.g. Maple, Mathematica). Moreover, it does not represent a good approach for languages consisting of words of fixed size (e.g. representing the search space of a combinatorial problem for a given instance), since the associated generating functions typically have very complex close-form expressions (if any), are hard to manipulate, and have infinite radius of convergence. In such cases, numerical iterations represent the method of choice.
<table><caption>Table 4.2: Weighted generators for context-free constructors, with $\pi(t)$ denoting the weight associated with a terminal letter/atom $t$. Remark that, in the case of products, the convolution can be restricted to values of $i$ in $[1, n-1]$ if the grammar is in $\varepsilon$-free form. Such a form can be assumed without loss of generality, and greatly eases the implementation by avoiding ill-defined equations, associated with non-productive circular dependencies in the precomputation.</caption><thead><tr><th>Rule type</th><th>Precomputation</th><th>Recursive generator ΓV<sub>n</sub></th></tr></thead><tbody><tr><td>Product v → a.b</td><td>v<sub>n</sub> ← Σ<sup>n</sup><sub>i=0</sub> a<sub>i</sub> × b<sub>n-i</sub></td><td>Choose K with prob.<br>P(K = k) = <sup>a<sub>k</sub> × b<sub>n-k</sub></sup>/<sub>v<sub>n</sub></sub><br><b>Return</b> ΓA<sub>k</sub>.ΓB<sub>n-k</sub></td></tr><tr><td>Union v → a, v → b</td><td>v<sub>n</sub> ← a<sub>n</sub> + b<sub>n</sub></td><td><b>If</b> Bernoulli(a<sub>n</sub>/v<sub>n</sub>) : <b>Return</b> ΓA<br><b>Else</b> : <b>Return</b> ΓB<sub>x</sub></td></tr><tr><td>Terminal v → t</td><td>v<sub>n</sub> ← <span style="font-size: 2em; vertical-align: middle;">{</span> <table style="display: inline-table; vertical-align: middle; border: none; text-align: center;"><tbody><tr><td>π(t) if n = 1</td></tr><tr><td>0 otherwise</td></tr></tbody></table></td><td><b>Return</b> t</td></tr><tr><td>Epsilon v → ε</td><td>v<sub>n</sub> ← <span style="font-size: 2em; vertical-align: middle;">{</span> <table style="display: inline-table; vertical-align: middle; border: none; text-align: center;"><tbody><tr><td>1 if n = 0</td></tr><tr><td>0 otherwise</td></tr></tbody></table></td><td><b>Return</b> ε</td></tr></tbody></table>

**Weighted random generation.** Generating within a weighted distribution is arguably the easiest part of multidimensional Boltzmann sampling.

Using the **recursive method**, weights can straightforwardly incorporated in classic recurrences, allowing to precompute the weights accessible through local alternatives. Those weights can, in turn, be used to derive probabilities that induce the weighted distribution described in Equation (4.9). Table 4.2 summarizes some of the constructors used in the context of context-free languages. The computation of those recurrences can be optimized using properties of holonomic generating functions, leading to an overall complexity in $O(mn \log n + n^{1+o(1)})$ for the generation of m words of length n in the weighted distribution of Equation (4.9), as discussed in Denise et al [50].

**Weighted Boltzmann samplers** can also be adapted from their uniform counterparts. To that purpose, it suffices to substitute, in the generators of Table 4.1, **weighted generating functions** defined as

$$C_{\pi}(z) = \sum_{\omega \in \mathcal{C}} z^{|\omega|} \prod_{i=1}^{k} \pi_{i}^{|\omega|_{t_{i}}}, \qquad (4.11)$$

instead of the ordinary ones used within classic Boltzmann sampling. As shown in Denise et al [50], such weighted generating functions are the solutions of systems of functional equations that is almost identical to those of Table 4.1, except for the addition of weights within **terminal rules** $v \to t$, namely

$$(\text{Ordinary gen. fun.}) C(z) = z \quad \to \quad (\text{Weighted gen. fun.}) C_{\pi}(z) = \pi(t) \times z \quad (4.12)$$

where $\pi(t)$ denotes the weight of an atom $t$. For any given weight vector, the Boltzmann parameter calibration can
### 4.2.2.3 Average-case complexity analysis

A usual with rejection-based samplers, the worst-case complexity is infinite, so we investigate the average-case complexity as the relevant metric to assess the efficiency of a sampler. In general, the probability to draw an object $\omega \in C$ having suitable composition

$$ (|\omega|_{t_1} = n \cdot r_1, |\omega|_{t_2} = n \cdot r_2 \dots) \in \mathbb{N}^k $$

in the atoms $(t_1, t_2 \dots)$ is given by

$$ \mathbb{P}(|W|_{t_i} = n \cdot r_i, \forall i \in [1, k] \mid |W| = n) = \frac{\sum_{\substack{\omega \in C \\ \text{s.t. } |\omega|=n, \\ |\omega|_{t_i} = r_i \cdot n, \forall i \in [1, k]}} \prod_{i=1}^{k} \pi_i^{|\omega|_{t_i}}}{\sum_{\substack{\omega|=n \\ \text{s.t. } |\omega|=n}} \prod_{i=1}^{k} \pi_i^{|\omega|_{t_i}}} \quad (4.13) $$

$$ = \frac{[z^n u_1^{n \cdot r_1} u_2^{n \cdot r_2} \dots] C_\pi(z; \pi_1 \cdot u_1, \pi_2 \cdot u_2 \dots)}{[z^n] C_\pi(z; \pi_1, \pi_2 \dots)} \quad (4.14) $$

where $C_\pi(z; u_1, u_2 \dots)$ represents the multivariate generation function introduced in Equation (4.10). The expected number of attempts is simply obtained as the inverse of the probability of acceptance, i.e. by exchanging the numerator and denominator in Equation (4.14).

**Simple type grammars.** For simple type grammars, associated with strongly connected aperiodic specifications, the multivariate version of the Drmota Theorem [56] applies to weighted specifications [50], and the composition in the k elected atoms follows a k-dimensional Gaussian limiting distributions, centered on the expectation of Equation (4.7) [56]. More precisely, the coefficients of the weighted multivariate generating function admit the following asymptotic equivalent:

$$ [z^n u_1^{m_1} u_2^{m_2} \dots] C_\pi(z; \pi_1 u_1, \pi_2 u_2 \dots) = \frac{\kappa_{\pi_1, \pi_2} \dots \rho_{\pi_1, \pi_2}^{-n}}{2^{k/2} n^{k/2} n \sqrt{n}} \left( e^{-n \Xi_{\pi_1, \pi_2} \dots} + O\left(\frac{1}{\sqrt{n}}\right) \right) $$

where $\kappa_{\pi_1, \pi_2} \dots$ is a constant of both n and k, $\rho_{\pi_1, \pi_2} \dots$ is the dominant singularity, and

$$ \Xi_{\pi_1, \pi_2} := \sum_{i=1}^{k} \sum_{j=1}^{k} \delta_{\pi_1, \pi_2}^{(i,j)} \times (\mu_i - m_i/n) \times (\mu_j - m_j/n). \quad (4.15) $$

The $\Xi_{\pi_1, \pi_2} \dots$ term essentially grows quadratically with k, and depends both on the deviation of the requested composition to the expected one $(\mu_1, \mu_2 \dots)$, and of values $\delta_{\pi_1, \pi_2}^{(i,j)} \dots$ that are related to the variance/co-variance matrix of the system. Remark that, upon a successful clibration of the weights $(\mu_i = r_i)$, $\Xi_{\pi_1, \pi_2} \dots$ vanishes at $(m_i = n \cdot r_i)_{i=1}^k$ so one has

$$ [z^n u_1^{n \cdot r_1} u_2^{n \cdot r_2} \dots] C_\pi(z; \pi_1 u_1, \pi_2 u_2 \dots) = \frac{\kappa_{\pi_1, \pi_2} \dots \rho_{\pi_1, \pi_2}^{-n}}{2^{k/2} n^{k/2} n \sqrt{n}} \left( 1 + O\left(\frac{1}{\sqrt{n}}\right) \right) $$
Moreover from the **Drmota-Lalley-Woods Theorem** [56, 112, 202] (see Theorem ??) we know that similar conditions lead a **square root type** singularity to dominate the ordinary generating function, inducing the following *universal* asymptotic equivalent for the coefficients

$$[z^n] C_\pi(z; \pi_1, \pi_2, \dots) = \frac{\lambda_{\pi_1, \pi_2, \dots} \rho_{\pi_1, \pi_2, \dots}^{-n}}{n\sqrt{n}} \left(1 + O\left(\frac{1}{\sqrt{n}}\right)\right).$$

The expected number of attempts probability of drawing an object of suitable composition, assuming successful calibration is then asymptotically equivalent to

$$\mathbb{E}_{\pi_1, \pi_2, \dots}(T) = \frac{[z^n] C_\pi(z; \pi_1, \pi_2, \dots)}{[z^n u_1^{n r_1} u_2^{n r_2} \dots] C_\pi(z; \pi_1 u_1, \pi_2 u_2, \dots)} \quad (4.16)$$

$$= 2^{k/2} n^{k/2} \left( \frac{\lambda_{\pi_1, \pi_2, \dots}}{\kappa_{\pi_1, \pi_2, \dots}} + O\left(\frac{1}{\sqrt{n}}\right) \right) \in \Theta\left(2^{k/2} n^{k/2}\right). \quad (4.17)$$

### 4.2.3 Beyond simplicity

The validity of those complexities extends beyond simple type grammars, and applies to any specification whose associated dominant singularity is unique, and arises from the contribution of a single connected component for the chosen weight assignment. It also applies to **periodic systems** such that, beyond a certain size N

$$[z^n] C(z) \neq 0 \Leftrightarrow n \equiv 0 [p] \text{ for all } n > N$$

for some period p > 1 through a simple rewriting of the grammar, coupled with a suitable change of variables.

However, general extensions of those complexities beyond strongly-connected systems appear very challenging. Indeed relaxing the connectivity constraint gives rise to many counter-examples, and may even induce bimodal distributions that cause an exponential dependency on n for the number of attempts. In fact, we were able to design an algorithm which constructs and return a specification whose limiting distribution for a dedicated atom follows an arbitrary discrete (univariate) distribution, specified by a user as part of the input [12].

**Remark 4.2.2:** Empirically, in the worst case scenario of an abysmally small probability of acceptance, estimates for the expected waiting time can be refined during the execution of the rejection algorithm. Upon reaching the conclusion of impending doom, an exact recursive algorithm [49] can then be elected.
## 4.3 Applications of multidimensional sampling

As noted above, failure to satisfy simplicity only (potentially) affects the complexity, but does not bear any consequence on correctness (uniformity). One can then take a *leap of faith*, and adopt the multidimensional Boltzmann paradigm in contexts where the concentration of the distribution is not guaranteed in general (but typically does hold).

In the context of RNA design (see Chapter 3 for further details), we argue for a change of paradigm, where design is approached as a random generation problem [88, 90, 117, 158, 160, 210]. Our methods generally rely on a polynomial-time dynamic programming algorithm to compute the **dual partition function** of a target structure S*

$$Z'_{S^*} = \sum_{\substack{\text{Sequence w such that} \\ \text{w compatible with S^*}}} e^{-\beta \cdot E_{w,S^*}}$$

where $E_{w,S^*}$ represents the free-energy of S* assuming a nucleotide sequence w, and $\beta$ represents an arbitrary constant. A **stochastic backtrack** algorithm reinterprets the recurrences into a rewriting system allowing the random generation of a compatible with emission probability

$$P(w | S^*) = \frac{e^{-\beta \cdot E_{w,S^*}}}{Z'_{S^*}}.$$

### 4.3.1 Controlling the GC-content (88, 90, 158, 160, 192, 193)

The **GC-content** of a sequence w simply represents its proportion of Guanines and Cytosines, defined as

$$GC(w) = \frac{|w|_G + |w|_C}{|w|}.$$

It is observed to vary substantially (16% to 75%) across genomes [94], or even between coding and non-coding regions within certain species [99], resulting from a selective pressure on the overall stability of RNAs.

Earlier methods for RNA design did not identify this parameter as an objective, leading to a strong bias towards high GC-contents observed among designed sequences, as revealed by our studies [158, 160]. Even more troublesome, was the observation that the GC-content of designed sequences typically appears Normal-distributed with linear variance. This means that a rejection strategy from the output of those tools would yield exponential runtime whenever the targeted GC-content does not coincide with the average GC-content of their output.

Fortunately, a multidimensional sampling strategy can be used to *shift* the expected GC-content of designed sequences to a given target. To that purpose, the dynamic program-
ming recurrences used for computing $Z'$ can be adapted through a minor modification to compute

$$Z'_{\gamma,S^*} = \sum_{\substack{w \\ \text{Sequence } w \text{ such that} \\ w \text{ compatible with } S^*}} e^{-\beta \cdot E_{w,S^*}} \times \gamma^{\text{GC}(w)}$$

where $\gamma$ is simply a positive weight (akin to $\pi_t$ above). A mirrored change in the stochastic backtrack allows to sample the GC-weighted distribution

$$P_{\gamma}(w | S^*) = \frac{e^{-\beta \cdot E_{w,S^*}} \times \gamma^{\text{GC}(w)}}{Z'_{\gamma,S^*}}.$$

This allows to shift the expected GC-content to a targeted GC-content $g \in [0, 1]$ by first calibrating $\gamma$, e.g. using a binary search, and then perform rejection sampling to restrict the output to sequences having a targeted GC-content.

The complexity analysis here becomes even harder to define than the above average-case analysis, as the GC distribution now depends on the target structure $S^*$. One way to measure the complexity is to consider, for any size $n$, the **worst average-case complexity** over all possible structures over $n$ nucleotides. In this setting, we could not formally establish the asymptotic complexity, but strongly suspect the expected number of attempts to grow like $O(\sqrt{n})$ in the worst-case. Indeed, popular energy models assign free-energy contributions to **loops** depending on their content, meaning that the GC contributions of two non-adjacent loops can be treated as independent random variables. We thus expect such limited, non-transitive, dependencies to allow an application of the **Central Limit Theorem**, and conclude on the normal distribution of the GC-content (sum of $\Theta(n)$ almost-independent variables).

Similar strategies were used for controlling the GC-content in the design of algorithms for probing evolutionary landscapes [192, 193], and for advanced versions of design requiring compatibility with multiple targets [88, 90].

### 4.3.2 Dinucleotides content of protein-coding RNAs (209).

Since the potential of sequences to adopt a given thermodynamic behavior cannot be easily captured at the generation stage, background models for RNAs usually include the conservation of a given **dinucleotides (2-mers) content**, defined for a sequence $w$ as the vector

$$d(w) = (d_{AA}(w), d_{AC}(w), d_{AG}(w), d_{AU}(w), \dots, d_{GG}(w))$$

where $d_{xy}(w)$ denotes the number of occurrences in $w$ for the 2-mer $xy \in \{A, C, G, U\}^2$. This property was essentially elected as an invariant, as it was empirically shown to act as a good proxy for the thermodynamic stability of RNAs [37]. Random RNA sequences
preserving a given 2-mers content can be uniformly generated as random Eulerian path, using the **dinucleotide shuffling algorithm** of Altschul and Erickson [6].

However, in the context of studying a protein-coding RNA m, a background model should ideally capture all known properties, to assess the statistical significance of a given phenomenon. Therefore, in addition to the 2-mers content, emitted sequences should also preserve the **reference protein** $P(m)$, encoded as an **amino-acids sequence**. This constraints allows for substitutions of **codons** (non-overlapping 3-mers), each encoding an amino-acid, such that the resulting sequence m' is said to be **synonymous** to m. In short, given a coding sequence m, $|m| \equiv 0$ [3], we must generate, uniformly at random, one or several RNA sequences synonymous to m' such as

$$d(m) = d(m') \quad \text{and} \quad P(m) = P(m').$$

Here, $P(\cdot)$ denotes the conversion of a coding RNA into a sequence of amino acids, achieved by a simple morphism

$$P(xyz.w) \to AA_{xyz}.P(w) \quad \text{and} \quad P(\varepsilon) \to \varepsilon$$

where $AA_{xyz}$ represents the amino-acid associated with the codon xyz $\in \{A, C, G, U\}^3$.

Due to the highly non-sequential nature of the Altschul and Erickson algorithm [6], there is little hope of adapting it to enforce synonymous substitutions. However, synonymous sequences can be modeled as a simple linear language $\mathcal{L}(m)$, defined by

$$\mathcal{L}(\varepsilon) = \varepsilon \quad \text{and} \quad \mathcal{L}(xyz.w) = \{x'y'z' | AA_{x'y'z'} = AA_{xyz}\} \times \mathcal{L}(w).$$

A predefined 2-mers content can then be enforced using multidimensional Boltzmann sampling, using 15 dimensions to represent the number of occurrences of each but one of the 2-mers (the last one being indirectly captured as the complement to the length). A weights vector $\mu := (\mu_{xy})_{xy \in \{A, C, G, U\}^2}$ is then introduced, defining a weighted distribution

$$\mathbb{P}_{\mu}(m' | m) = \frac{\prod_{xy \in \{A, C, G, U\}^2} \mu_{xy}^{|m'|_{xy}}}{\mathbb{Z}_{\mu, m}} \quad (4.18)$$

where

$$\mathbb{Z}_{\mu, m} \equiv \sum_{\substack{m' \text{ s.t.} \\ P(m') = P(m)}} \prod_{xy \in \{A, C, G, U\}^2} \mu_{xy}^{|m'|_{xy}}$$

The above partition function can be computed in linear-time using dynamic programming

$$\mathbb{Z}_{\mu,abc.w;\alpha} := \sum_{\substack{\alpha'b'c' \in \{A,C,G,U\}^3 \\ \text{s.t. } AA_{\alpha'b'c'} = AA_{abc}}} \mathbb{Z}_{\mu,w;c'} \times \prod_{xy \in \{A,C,G,U\}^2} \mu_{xy}^{\alpha abc|_{xy}}$$

with $\mathbb{Z}_{\mu,\varepsilon;\alpha} := 1$ and $\mathbb{Z}_{\mu,m} \equiv \mathbb{Z}_{\mu,m;\varepsilon}$. A stochastic sampling algorithm can then be adapted in a straightforward manner to sample with respect to the distribution (4.18),
and coupled with a rejection scheme to obtain the algorithm at the core of the SPARCS method [209]. Complexity-wise, the 15 degrees of freedom induce an expected $\Theta(n^{15/2})$ attempts for each admissible sequence, under the assumption of a Normal multidimensional distribution. In practice, we found that a tolerance of a couple of occurrences for each 2-mers was sufficient to induce a linear empirical runtime for the method.

## 4.4 Redundancy in weighted sampling and countermeasures

Redundancy within of sampled sets of objects constitutes a possible shortcoming of methods relying on a weighted sampling. Indeed, as argued in Section 2.3.3, redundancy is generally useless in inferential statistics whenever the emission probability of a sample can be exactly determined a posteriori. However, the extent to which redundancy prevails quantitatively within random samples, and thus has a measurable adverse effect on the efficiency of estimators, had not been the object of precise inquiries beyond the seminal contribution of Flajolet, Gardy and Thimonier [70].

Let us first define formally the key concept of redundancy.

**Definition 4.4.1 (Redundancy of sample):** Given a collection $\Omega$, the redundancy of a sample $(e_1, \dots, e_M) \in \Omega^M$ is defined as

$$ \Re(e_1, e_2, \dots, e_M) = \frac{M}{|(e_1, e_2, \dots, e_M)|} - 1 \quad (4.19) $$

### 4.4.1 Collisions in weighted languages (58, 78)

The study of the expected redundancy is very much connected to the **waiting times** of several events occurring in probabilistic allocations. The **Birthday Paradox** (BP) problem represents one such instance, and considers the time of the first collision.

**Definition 4.4.2 (Birthday paradox):** Given a collection $\Omega$ equipped with a probability distribution, the **Birthday** $B_\Omega$ is defined as

$$ B_{\Omega} = \mathbb{E}(M | \Re(E_1, \dots, E_{M-1}) = 0 \wedge \Re(E_1, \dots, E_M) > 0) $$

where $E_i$ are random independent variables in $\Omega$.
For instance, in the uniform distribution over $\Omega$, one has

$$ \mathbb{P}(\text{Collision at time } i \mid \{e_1, \dots, e_{i-1}\}) = \frac{|\{e_j\}_{j=1}^{i-1}|}{|\Omega|} < \frac{B_{\Omega}}{|\Omega|}, \quad \forall i \le B_{\Omega} $$

and it follows that, after $B_{\Omega}$ generations, one has

$$ \mathbb{E}(\#\text{Collisions}) < \frac{B_{\Omega} \times B_{\Omega}}{|\Omega|} \in \mathcal{O}(1) $$

and we conclude that $\mathbb{E}(\Re(e_1, \dots, e_{B_{\Omega}})) \to 0$ for asymptotically large collections $\Omega$.

Another setting of interest is that of the **Coupon Collector** (CC) problem, the time which considers the expected sample size before the full collection is obtained.

**Definition 4.4.3 (Coupon Collector):** Given a collection $\Omega$ equipped with a probability distribution, the **Coupon Collector** $C_{\Omega}$ is defined as

$$ C_{\Omega} = \mathbb{E}(M \mid \{e_1, \dots, e_M\} = \Omega) \quad (4.20) $$

It can be shown that the expected redundancy is an increasing function of $M$, so it is natural to bound the redundancy is by that of the full collection, namely

$$ \mathbb{E}(\Re(e_1, \dots, e_M) \mid \{e_1, \dots, e_M\} = \Omega) = \frac{C_{\Omega}}{|\Omega|} - 1. $$

In collaboration with Danièle Gardy and Jérémie du Boisberranger [58, 78], we revisited these quantities attempting to obtain precise specialization of the general results of Flajolet, Gardy and Thimonier [70] in the specific context of **weighted languages**.

### 4.4.1.1 Main results

Let $C$ be a specification generating a language $\mathcal{L}$, with a collection $\pi$ of weights associated with atoms. We denote by $\pi_n^{\nabla}$ (resp. $\pi_n^{\nabla}$) the **minimal** (resp. **maximal**) **weight of a word** in $\mathcal{L}_n$. Consider the weighted generating function

$$ C_{\pi}(z) = \sum_{\omega \in \mathcal{L}} \pi_{i} z^{|\omega|} = \sum_{n \ge 0} \mu_{\pi, n} z^{|\omega|} $$

and denote by $\rho_{\pi}$ its dominant singularity. Finally let us define the $k$-th **moment** of a $\pi$-weighted distribution as

$$ \alpha_{k,n} = \sum_{i=1}^{m_n} p_i^k = \frac{\sum_{w \in \mathcal{L}_n} \pi(w)^k}{\mu_{\pi,n}^k} = \frac{\mu_{\pi^k, n}}{\mu_{\pi, n}^k}. \quad (4.21) $$

We consider the following conditions:
**C1 Diversity:** The max. probability $\pi_n^\Delta/\mu_{\pi,n}$ of a word in $\mathcal{L}_n$ is exponentially decreasing with n;

**C2 Log-positive weights:** For each terminal symbol $t \in \Sigma$, $\pi_t > 1$;

**C3 Bounded dependency:** For any rational number $k > 1$ and any weight vector $\pi$ such that Condition C2 holds, $\rho_\pi^k < \rho_{\pi^k}$ holds.

**Theorem 6 (Birthday paradoxes in weighted languages):** Under conditions C1, C2 and C3, the birthday paradox of $\mathcal{L}_n$ is given by

$$ \mathcal{B}_{\mathcal{L}_n; \pi} \sim \frac{\sqrt{\pi}}{\sqrt{2\alpha_{2,n}}} = \frac{\mu_{\pi,n}\sqrt{\pi}}{\sqrt{2\mu_{\pi^2,n}}} \in \Omega(\gamma^n), \quad \gamma := \frac{\rho_{\pi}}{\sqrt{\rho_{\pi^2}}} > 1 $$

In other words, the first collision typically occurs after a time exponential on n (albeit with typically small base), inducing essentially zero redundancy before $\mathcal{B}_{\Omega;\pi}$ samples.

**Theorem 7 (Coupon collector for weighted languages):** The expected number of generations $\mathcal{C}_{\mathcal{L}_n;\pi}$ needed to get all words in $\mathcal{L}_n$ is such that

$$ \frac{\mu_{\pi,n}}{\pi_n^\nabla} \le \mathcal{C}_{\mathcal{L}_n;\pi} \le 2 \cdot \mathcal{H}_{m_n} \cdot \frac{\mu_{\pi,n}}{\pi_n^\nabla} \quad (4.22) $$

This result is a direct corollary of Berenbrink and Sauerwald [16], and was later refined in Du Boisberranger et al. [58]. In particular, Equation (4.22) implies that, for a vast collection of languages the **expected redundancy** can grow exponentially with n for a given targeted occupancy.

**Remark 4.4.1 (Redundancy of RNA sampling (homopolymer model)):** In particular, this results impacts any specifiable family of Boltzmann-weighted RNA secondary structures (homopolymer models), for which one typically has $\pi_n^\nabla \in \mathcal{O}(1)$. It implies that, when the full collection is reached, the **expected number of copies** of a secondary structure essentially grows (up to a possible $\Theta(n)$ factor) like $\mu_{\pi,n}/|\mathcal{L}_n| \in \Theta(\alpha^n)$ with $\alpha > 1$. This forbids using rejection-based generation strategies for the non-redundant generation of large sample sizes, since rejections would then become overwhelming.

Moreover, denote by $W_{n,\pi}$ the **weight classes** in $\mathcal{L}_n$, such that

$$ W_{n,\pi} := \{(w, m) \mid w \in \mathcal{L}_n : \pi(w) = w \text{ and } m = \{|w \in \mathcal{L}_n | \pi(w) = w|\} \} $$

Then the expected **number of distinct elements after k generations** can be computed.
**Theorem 8 (Distinct weighted words):** The **expected number of distinct words** after k generations is given by

$$ \mathbb{E}(|\{e_1, \dots, e_k\}|) = \sum_{(w,m) \in W_{n,\pi}} m \cdot \left(1 - \left(1 - \frac{w}{\mu_{\pi,n}}\right)^k\right) = \sum_{(w,m) \in W_{n,\pi}} m \cdot \left(1 - e^{-\frac{w}{\mu_{\pi,n}}k}\right) + \mathcal{O}(1). $$

A similar results holds for the **coverage C** of a sample $(e_1, \dots, e_k)$, defined as

$$ C(e_1, \dots, e_k) = \frac{\sum_{e \in \{e_1, \dots, e_k\}} \pi(e)}{\mu_{\pi,n}} \quad (4.23) $$

which represents the proportion of the distribution that has already been sampled.

**Theorem 9 (Coverage of weighted words):** The **expected coverage** of a set of distinct words after k generations is given by

$$ \mathbb{E}(C(e_1, \dots, e_k)) = \sum_{(w,m) \in W_{n,\pi}} m \cdot \frac{w}{\mu_{\pi,n}} \cdot \left(1 - \left(1 - \frac{w}{\mu_{\pi,n}}\right)^k\right). $$

Since there are at most $(n+1)^{|\Sigma|}$ different compositions/classes of weights, the expressions of Theorems 8 and 9 can be computed in polynomial time algorithms from the **weighted-classified partition functions** $\mu_{\pi,n}$, also known as the **density of states**.

Since redundancy can become overwhelming, it is natural to pursue algorithmic countermeasures. More precisely, the above result shows that the algorithmic problem of sampling without replacement from a weighted language cannot be efficiently tackled using rejection, motivating an unranking approach [118], or a modified instance of the recursive method [118, 147], illustrated in Section 2.3.3 and used to build RNA kinetics landscapes [133].
Chapter 5

Asymptotic properties of RNA secondary
structures and other trees

The tree-like nature of RNA secondary structures not only makes it an excellent play-
ground for dynamic programming, but also enables the study of its generic asymptotic
properties. Such work go back to the very early days of RNA bioinformatics, with pioneer-
ing works from Waterman [194], Zuker and Sankoff [212], Viennot and Vauchaussade de
Chaumont [187], Hofacker [96], Bundschuh and Hwa [24], Nebel [139], Clote [36], Reidys
and Jin [104], and others too numerous to mention exhaustively. While initially moti-
vated by algorithmic [194] and evolutionary [212] considerations, the field has gradually
turned towards questions in statistical mechanics, addressing the behavior of RNA at the
thermodynamic equilibrium, as exemplified by the work of Bundschuh and Hwa [25]. In
this context, asymptotic enumerations can be used to elucidate the equilibrium properties
of random large homopolymers, including their dependency to changes of temperatures
(phase transitions [23, 26]).

In the **homopolymer model**, nucleotides are anonymized, so that any pair can form **base pairs**, inducing structural motifs that are equally contributing to the stabilization of their conformation. Implicitly, this amounts to coarsely abstract an RNA as its length n. While such a simplification may appear bold at first sight, many of the properties derived from this model are robust to the incorporation/perturbations of the energy model, or can be used as a first step to recover the dependency on an actual nucleotide content (heteropolymer model) [24, 212].

This chapter summarizes some contributed asymptotic analyses, using techniques in
enumerative and analytic combinatorics, to address questions in algorithmic analysis
and design [33, 148, 150], RNA thermodynamics [38], network properties [180], and
design/evolution [206]. In the interest of self-completeness, we start by reminding basic
tools and notions. A magisterial exposition, including much more versatile techniques,
can however (and should!) be sought in the Flajolet/Sedgewick *bible* [69].
**Outline.** In this short chapter, I first remind in Section 5.1 some basic concepts, tools and principles in analytics combinatorics and, in Section 5.2, I describe the application of those techniques to selected problems involving the secondary structure of RNA.

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

W. A. Lorenz, P. Clote, and Y. Ponty. Asymptotics of RNA shapes. **Journal of Computational Biology**, 15(1): 31–63, 2008

Y. Ponty. Efficient sampling of RNA secondary structures from the Boltzmann ensemble of low-energy: The boustrophedon method. **Journal of Mathematical Biology**, 56(1-2):107–127, 2008

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

P. Clote, Y. Ponty, and J.-M. Steyaert. Expected distance between terminal nucleotides of RNA secondary structures. **Journal of Mathematical Biology**, 65(3):581–99, Sept. 2012

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

D. Surujon, Y. Ponty, and P. Clote. Small-world networks and RNA secondary structures. **Journal of computational biology: a journal of computational molecular cell biology**, 26(1):16–26, Jan. 2019

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

# 5.1 Basic tools in enumerative and analytic combinatorics

Analyzing the asymptotic enumerative behavior of a combinatorial class can be done by the following four-step program, referred to as the **symbolic method** in the pioneering work of the late Philippe Flajolet:

1. Find a unambiguously grammar/specification that generates the objects;

2. Translate the grammar into a system of functional equations;

3. Solve the system and identify the dominant singularities;

4. Use some automatic Theorem to extract the leading term for the asymptotic equivalent for the number of objects of size n.

## 5.1.1 Context-free languages.

A **combinatorial class** is a, possibly infinite, set $\mathcal{A}$ where every object $\alpha \in \mathcal{A}$ has a **length** $|\alpha|$, and such that the restriction $\mathcal{A}_n$ to objects of $\mathcal{A}$ having any given length $n$ is finite.
Combinatorial classes are usually described by a specification, which can be interpreted as a generative process. Among such processes are **context-free grammars**, abstract derivation systems formally defined as shown below.

**Definition 5.1.1 (Context-free grammar):** A Context-free grammar is a 4-tuple $G = (\mathcal{N}, \Sigma, \mathcal{R}, S_0)$ where:

* $\mathcal{N}$ is a set of intermediate symbols called **non-terminal symbols**;

* $\Sigma$ is a set of letters, also called **terminal symbols**, used for encoding our combinatorial objects;

* $\mathcal{R}$ is a set of **productions rules** of the general form $S \to w$, where $v \in \mathcal{N}$ is a non-terminal symbol, and $w$ is a word composed of non-terminal or terminal symbols. When the right-hand side is empty ($|w| = 0$), the rule is said to generate the **empty word** $\varepsilon$;

* $S_0$ is the **axiom**, a special non-terminal from which the generation is initiated.

A **derivation** results from the application of a production rule $S \to w$ in $\mathcal{R}$, replacing the non-terminal $S$ on the left-hand side with the right-hand side $w$. Starting from $w := S_0$, one can iterative choosing a production $S' \to w'$ for the leftmost non terminal $S'$ in $w$, replacing its occurrence with $w'$, until $w$ only consists of terminal letters only. One then obtains a tree-like sequence of choices, called a **parse tree**. Considering (implicitly) all possible parse trees leads to a, possibly infinite, set $\mathcal{L}_G$ of sequences, which is called the **language** of $G$ and can be interpreted as a combinatorial class.

**Definition 5.1.2 (Generated language):** The language $\mathcal{L}_S$ of a non-terminal $S \in \mathcal{N}$ within a context-free grammar $G = (\mathcal{N}, \Sigma, \mathcal{R}, S_0)$ is the set $\mathcal{L}_S$ such that

$$ \forall t \in \Sigma : \mathcal{L}_t = \{t\} \quad \text{and} \quad \forall S \in \mathcal{R} : \mathcal{L}_S = \bigcup_{(S \to w) \in \mathcal{R}} \prod_{i=1}^{|w|} \mathcal{L}_{w_i} \quad (5.1) $$

where $\prod$ is the Cartesian product of languages, i.e. all ordered concatenations of words from the sets, defined as $\{\varepsilon\}$ for an empty word $w = \varepsilon$.

The **language** of a CFG grammar $G$ is defined as $\mathcal{L}_G := \mathcal{L}_{S_0}$.

Any grammar can be transformed into an equivalent **Chomsky Normal Form** (CNF) grammar, where the rules associated to each non-terminal $S \in \mathcal{N}$ are restricted to the following cases:

* **Product type** non-terminal $S \to T.U$ where $T, U \in \mathcal{N}$;

* **Union type** non-terminal $S \to T$ and $S \to U$, where $T, U \in \mathcal{N}$;
* Terminal type non-terminal $S \to t$, where $t \in \Sigma$ is a terminal symbol;

* Epsilon type non-terminal $S \to \epsilon$, denoting the empty word;

Next comes the crucial notion of **unambiguity** of a grammar, which guarantees a unique way to produce each word in $L_G$.

**Definition 5.1.3 (CFG Unambiguity):** CFG $G = (\mathcal{N}, \Sigma, \mathcal{R}, S_0)$ is **unambiguous** if and only if one of the following (equivalent) statements is true:

* there exists a bijection between the parse trees and the elements of $L$;

* all parse trees produce distinct words over the terminal alphabet;

* the unions in Equation 5.1 are disjoint, and any word produced by a Cartesian products can be uniquely factored into its parts.

Another way if For the rest of this exposition, we are going to assume that the considered grammars are unambiguous and. Moreover, this time without loss of generality, we will assume that our grammars do not have **infinite cycles of unit productions**:

$$S \to T \to \dots \to S$$

Those are non productive and can be detected and eliminated by standard algorithms

## 5.1.2 Enumeration and generating functions.

Within unambiguous grammars, one can easily characterize the set of words $L_S$ generated from each non-terminal $S \in \mathcal{N}$ and the number $s_n$ of its elements of length $n$, depending of its type:

* Product $S \to T.U$: Each word generated from $S$ can be decomposed into elements in $L_T$ and $L_U$ of length $i$ and $n - i$ respectively. The grammar is unambiguous, so this **decomposition is unique**, so we have

$$s_n = \sum_{i=0}^{n} t_i \times u_{n-i};$$

* Union $S \to T$ and $S \to U$: The unambiguity of the grammar implies that the **union is disjoint**, so we have $s_n = t_n + u_n$;

* Terminal $S \to t$: Generates $\{t\}$, so we have $s_1 = 1$ and $s_n = 0, \forall n \neq 1$;

* Epsilon $S \to \epsilon$: Generates $\{\epsilon\}$, so we have $s_0 = 1$ and $s_n = 0, \forall n \neq 0$.

We can now introduce the key tool of enumerative combinatorics.
**Definition 5.1.4** ((Counting) generating function): A **generating function** $A(z)$ for a combinatorial $\mathcal{A}$ is a formal power series that summarizes the cardinality information over $\mathcal{A}$:

$$A(z) = \sum_{a \in \mathcal{A}} z^{|a|} = \sum_{n \ge 0} a_n z^n$$

where $a_n = |\mathcal{A}_n|$ denotes the number of objects of size n in $\mathcal{A}$.

In the above definition, $z$ is a formal complex variable, and can be thought as a medium of useful information, later retrieved using a Taylor expansion at $z = 0$. We use the notation $[z^n] A(z)$ to represent the n-th coefficient of $A(z)$, noting that $[z^n] A(z) = a_n$.

Consider the generating function $S(z)$ associated with (language of) a non-terminal $S \in \mathcal{N}$. By definition, one has

$$S(z) = \sum_{w \in \mathcal{L}_S} z^{|w|} = \sum_{n \ge 0} s_n z^n.$$

The recursive nature of the grammar then induces a fairly simple recursive scheme, and the generating function $S(z)$ of a non-terminal $S$ obeys a remarkably simple **system of functional equations**:

* Product $S \to T.U$: We have

$$S(z) = \sum_{n \ge 0} s_n z^n = \sum_{n \ge 0} \left( \sum_{i=0}^{n} t_i \times u_{n-i} \right) z^n = T(z) \times U(z) \quad (5.2)$$

* Union $S \to T$ and $S \to T$: We have

$$S(z) = \sum_{n \ge 0} s_n z^n = \sum_{n \ge 0} (t_n + u_n) z^n = \sum_{n \ge 0} t_n z^n + \sum_{n \ge 0} u_n z^n = T(z) + U(z).$$

* Terminal $S \to t$: $S(z) = \sum_{w \in \mathcal{L}_S} z^{|w|} = z^{|t|} = z$;

* Epsilon $v \to \varepsilon$: $S(z) = \sum_{w \in \mathcal{L}_S} z^{|w|} = z^{|\varepsilon|} = z^0 = 1$.

Any unambiguous context-free grammar in CNF can thus be transformed into a system of functional equation. This system is **algebraic**, as Gaussian elimination can be used to express it as

$$P(z, S(z)) = 0 \text{ where } P \text{ is a polynomial in } S(z) \text{ and } z, \quad (5.3)$$

and its solutions are called **algebraic generating functions** by analogy.

### Example 18: DNA strings

DNA strings can be expressed as a language $\{A, C, G, T\}^*$, generated by the context-free grammar

$$S \to A.S \mid C.S \mid G.S \mid T.S \mid \varepsilon$$
This grammar can be transformed in the CNF grammar for illustrative purposes

$S \to S' | E$

$a \to A$

$S' \to T.S$

$c \to C$

$T \to T' | T''$

$g \to G$

$T' \to a | c$

$t \to T$

$T'' \to g | t$

$E \to \epsilon$

The transforms described above, coupled with minor simplifications, lead to the following system

$$S(z) = 4z \times S(z) + 1 = \frac{1}{1 - 4z}$$

which can be solved to obtain

$$S(z) = \frac{1}{1 - 4z} = \sum_{n \ge 0} 4^n z^n.$$

**Remark 5.1.1:** In general, algebraic systems associated with context-free grammars are non-linear, so that conjugate solutions may exist. However, only one will admit positive integers as its coefficients.

### 5.1.3 Basic singularity analysis

Finding a closed-form expression for $s_n$, the coefficients of $S(z)$, is generally a highly involved task (although such a form always exists [11]). An easier alternative is to use **analytic combinatorics** to extract a simple asymptotic equivalent for $s_n$. Singularity analysis simplifies this problem by focusing on the **singularities of** $S(z)$, i.e. points of the complex plane where the generating function ceases to be analytic.

Of particular interest are the **dominant singularities**, i.e. singularities of smallest modulus that are the main contributors to the coefficients. From Pringsheim's Theorem [69, pp 240], we know that at least one of them lies on the real axis for non negative series. Moreover, only limited types of singularities may occur in algebraic generating functions. This leads to predictable behaviors for the coefficients, and allows to bypass complex analysis (almost) entirely thanks to *automatic* theorems.

<table><caption>**Table 5.1:** Generated languages, enunerations and ordinary (counting) generating functions for unambiguous Chomsky Normal Form grammars.</caption><thead><tr><th>Rule type</th><th>Language</th><th>Coefficients</th><th>Generating function</th></tr></thead><tbody><tr><td>Product S → T.U</td><td>L<sub>S</sub> = L<sub>T</sub> × L<sub>U</sub></td><td>s<sub>n</sub> = &sum;<sub>i=0</sub><sup>n</sup> t<sub>i</sub> × u<sub>n-i</sub></td><td>S(z) = T(z) × U(z)</td></tr><tr><td>Union S → T, S → U</td><td>L<sub>S</sub> = L<sub>T</sub> ∪ L<sub>U</sub></td><td>s<sub>n</sub> = t<sub>n</sub> + u<sub>n</sub></td><td>S(z) = T(z) + U(z)</td></tr><tr><td>Terminal S → t</td><td>L<sub>S</sub> = {t}</td><td>s<sub>n</sub> = <table><tr><td>1</td><td>if n = 1</td></tr><tr><td>0</td><td>otherwise</td></tr></table></td><td>S(z) = z</td></tr><tr><td>Epsilon S → ε</td><td>L<sub>S</sub> = {ε}</td><td>s<sub>n</sub> = <table><tr><td>1</td><td>if n = 0</td></tr><tr><td>0</td><td>otherwise</td></tr></table></td><td>S(z) = 1</td></tr></tbody></table>
**Theorem 10 (Asymptotic growth (69, pp 244))**: Let $\rho \in \mathbb{R}^+$ be the dominant singularity of an algebraic generating function $S(z)$, then:

$$ \lim_{n \to +\infty} \sqrt[n]{[z_n] S(z)} = 1/\rho. $$

In other words, the exponential part of the asymptotics is solely driven by the value of the dominant singularity. This property is sometimes referred to as the **first principle of Analytic Combinatorics**, and is sufficient to obtain crude estimates.

### Example 19: DNA strings – Asymptotics

We illustrate this principle on DNA strings, associated with a generating function $S(z) = \frac{1}{1-4z}$. Here, the (unique) dominant singularity is a pole (cancellation of the denominator) at $z = 1/4$.

Theorem 10 implies that, for some subexponential term $f(n)$, $(\lim_{n \to +\infty} \sqrt[n]{f(n)} = 1)$,

$$ [z^n] S(z) = 4^n \cdot f(n). $$

An application of Theorem 11 below refines this estimate to $[z^n] S(z) = 4^n(1 + O(1/n))$, meaning that we can rule out the existence of an additional polynomial term.

The **second principle of Analytic Combinatorics** considers the *nature* of singularities to determine more precise estimates. **Transfer theorems** can then be used to relate the asymptotic growth of $s_n$ to the behavior of $S(z)$ close to its dominant singularity.

**Theorem 11 (Rational Singularities (69, pp 256))**: Let $\rho \in \mathbb{R}^+$ be the dominant singularity of an algebraic generating function $S(z)$, rewritten as

$$ S(z) = \frac{U(z)}{(1 - z/\rho)^\alpha} + T(z) \qquad (5.4) $$

where $\alpha \in \mathbb{N}^+ - \{0\}$, $\kappa \in \mathbb{R}$, $(1 - z/\rho)^\alpha \times T(z)$ and $U(z)$ are analytic at $z = \rho$ and $\kappa := U(\rho)$. Then the coefficients of $S(z)$ admit the following asymptotic equivalent

$$ [z^n] S(z) = \frac{\kappa}{(\alpha - 1)!} \times n^{\alpha - 1} \times \rho^{-n} (1 + O(1/n)). \qquad (5.5) $$

Polar singularities are the only singularities arising from left/right linear grammars (or, more generally, rational languages), associated with systems in Equation 5.3 that are linear in $S(z)$. In such cases, the solution of the system is a **rational generating function**, i.e. it can be reduced to a form $A(z)/B(z)$, for A and B some polynomials. The alternative formulation required by Equation (5.6), can thus always be found.

### Example 20: Fibonacci language
The Fibonacci language is generated by

$$F \to aF | bF | \varepsilon$$

The number $f_n$ of words of size $n$ is thus given by

$$f_n = \begin{cases} f_{n-1} + f_{n-2} & \text{if } n \ge 2 \\ 1 & \text{otherwise } (n \in \{0, 1\}) \end{cases}$$

such that $f_n$ is the $(n+1)$-th Fibonacci number. The generating function is then

$$F(z) = zF(z) + z^2 F(z) + 1 = \frac{1}{1-z-z^2} = \frac{1}{(1-z/\rho_+)(1-z/\rho_-)} \quad \text{with } \rho_\pm = \frac{-1 \pm \sqrt{5}}{2}.$$

Here, the dominant singularity is $\rho_+$ since $|\rho_+| < |\rho_-|$. $F(z)$ can then be reformulated as

$$F(z) = - \frac{\rho_- / \sqrt{5}}{1 - z / \rho_+} + \frac{\rho_+ / \sqrt{5}}{1 - z / \rho_-}$$

where the second term is already analytic at $z = \rho_+$ (and so is its multiplication by $1 - z/\rho_+$).

A direct application of Theorem 11, with $\rho := \rho_+$, $\alpha := 1$, and $\kappa := \rho_- / \sqrt{5}$, provides a precise asymptotic estimate

$$[z^n] F(z) = \frac{\rho_-}{\sqrt{5}} \rho_+^{-n} (1 + \mathcal{O}(1/n)).$$

Polar singularities may also dominate in (non linear) cases, and their simplicity always make them worthy of an initial investigation. However, non-rational languages will typically see more complex singularities arise and dominate their asymptotics. Fortunately, as covered by the following result, similar asymptotics hold, and can be derived at the cost of marginally more work.

**Theorem 12 (Single Algebraic Singularity (69, pp 393) + (11, 68))**: Let $\rho \in \mathbb{R}^+$ be the single dominant singularity of an algebraic generating function $S(z)$, rewritten as

$$S(z) = \frac{U(z)}{(1 - z/\rho)^\alpha} + T(z) \qquad (5.6)$$

where $\alpha \in \mathbb{Q} - \mathbb{N}$, $T(z)$ and $U(z)$ are analytic at $z = \rho$, and $\kappa := \lim_{z \to \rho} U(z) > 0$. Then the coefficients of $S(z)$ admit the following asymptotic equivalent

$$[z^n] S(z) = \frac{\kappa}{\Gamma(\alpha)} \times n^{\alpha-1} \times \rho^{-n} (1 + \mathcal{O}(1/n)).$$

The main difference with the rational case covered by Theorem 11, is that $\alpha$ can take more general values than just positive ones. This way, it covers singularities arising from radical forms $\sqrt{1 - z/\rho}$, called **square-root type singularities**, that are ubiquitous when enumerating tree-like objects.

The case of **multiple dominant singularities**, i.e. singularities $\rho_0, \rho_1, \rho_2, \dots$ such that $|\rho_0| = |\rho_1| = |\rho_2| = \dots$, is slightly more involved. In cases where the co-dominant
singularities lead to different values $\alpha_0 < \alpha_1 \le \alpha_2$, then only the largest needs to be considered, and Theorem 12 holds with a slightly enlarged error term $O(1/n^{\min(\alpha_0-\alpha_1,1)})$.

However, singularities are complex numbers and may team with, or cancel each others' contributions for certain values of $n$ (e.g. in the case of complex conjugates). However, as mentioned in Banderier and Drmota [11], multiple dominant singularities will give rise to periodic behaviors, with asymptotic equivalents still obeying the general form given in Theorem 12 (or 0) for any value of $n$ such that $n \equiv r[p]$, where $p$ is the **period** of the system.

In practice, a grammar featuring multiple dominant singularities can be rewritten, to only generate words of length $n$ with a prescribed modular value $r$ which, in combination with a suitable change of variables/vocabulary $\Sigma$, will suppress such periodic behaviors and lead to a single dominant singularity.

**Remark 5.1.2:** Banderier and Drmota [11] provide a definitive characterization of singularities in algebraic generating functions, including those of context-free languages. They show that the $\alpha$ value for dominant singularities can only take **dyadic values**

$$ \alpha \in \left\{ -\frac{1}{2^k} \mid k \ge 1 \right\} \cup \left\{ \frac{m}{2^k} \mid m \ge 1 \wedge k \ge 0 \right\} \subset \mathbb{Q}. $$

By further restricting the type of dominant singularities in algebraic generating functions, this results provides new tools to prove the intrinsic ambiguity of certain classes of languages and, more generally, the impossibility to model combinatorial classes using certain types of specifications.

## Example 21: RNA Secondary Structures

Following Waterman [194], RNA secondary structures can be generated by the grammar

$$ S \rightarrow (S^{\ge 1}) S \mid \bullet S \mid \varepsilon $$

where $S^{\ge 1}$ forbids $S$ to generate the empty word $\varepsilon$, and thus has generating function $S(z) - 1$. The associated system of functional equations is then

$$ S(z) = z(S(z) - 1)zS(z) + zS(z) + 1 \Leftrightarrow 0 = z^2S^2(z) + (-1 + z - z^2)S(z) + 1 $$

Solving the quadratic equation, gives two solutions for $S(z)$

$$ S^{\pm}(z) = \frac{1 - z + z^2 \pm \sqrt{\Delta}}{2z^2} $$

where $\Delta$ is the discriminant

$$ \begin{aligned} \Delta &:= (1 - z + z^2)^2 - 4z^2 = (1 + z + z^2)(1 - 3z + z^2) \\ &= \left(1 - \frac{z}{x_{-}}\right)\left(1 - \frac{z}{x_{+}}\right)\left(1 - \frac{z}{y_{-}}\right)\left(1 - \frac{z}{y_{+}}\right) \text{ with } x_{\pm} := \frac{3 \pm \sqrt{5}}{2} \text{ and } y_{\pm} := \frac{-1 \pm i\sqrt{3}}{2} \end{aligned} $$

Extracting the first terms of $S^+(z)$ reveals negative coefficients, inconsistent with our expectation of a non-negative number of structures, and we conclude that

$$ S(z) := S^{-}(z) = \frac{1 - z + z^2 - \sqrt{(1 - z/x_{-})(1 - z/x_{+})(1 - z/y_{-})(1 - z/y_{+})}}{2z^2}. $$
The function admits singularities for $z \in \{x_-, x_+, y_-, y_+\}$, and is analytic everywhere else in $\mathbb{C}$. The dominant singularity is $\rho := x_-$ as it has smallest modulus, and we can rewrite $S(z)$ as

$$S = \frac{1 - z + z^2}{2z^2} - \frac{1}{(1 - z/\rho)^{-1/2}} \times \mathcal{U}(z) \quad \text{with} \quad \mathcal{U}(z) := \frac{\sqrt{(1 - z/x_+) (1 - z/y_-) (1 - z/y_+)}}{2 z^2}$$

where $\mathcal{U}(z)$ can be verified to be analytic at $z = \rho$, and reaches the following value

$$\kappa := \mathcal{U}(\rho) = \frac{\sqrt{12\sqrt{5} - 20}}{(7 - 3\sqrt{5}) \times \sqrt{3} + \sqrt{5}} = \sqrt{\frac{15 + 7\sqrt{5}}{2}}.$$

A direct application of Theorem 12, with $\alpha = -1/2$, gives the final asymptotic equivalent:

$$[z^n] S(z) = \frac{\kappa}{2\sqrt{\pi}} \times \frac{\rho^{-n}}{n\sqrt{n}} (1 + \mathcal{O}(1/n))$$

This asymptotic first-order estimate is remarkably accurate, and is already within 1% of the true count for RNAs of length 250 nts, and keeps on getting better and better to reach asymptotic equality.

### 5.1.4 Useful extensions and shortcuts

#### 5.1.4.1 Ensemble properties using multivariate gen. fun. and derivatives

So far, generating functions have only provided valuable information on the number of words of size $n$ within a language. However, in many cases, the asymptotic growth of $s_n = [z^n] S(z)$ is arguably less relevant to a downstream application than the expected properties of a word in $\mathcal{L}_n$. To study such properties, one needs to (transiently) introduce a bivariate generating functions, which not only considers the length of an object, but also some additional property, here called a feature.

Namely, consider **atomic contributions** $\{f_r\}_{r \in \mathbb{R}}$, associating a (possibly null) numeric value to the production rules of a grammar.

**Definition 5.1.5 (Additive feature):** An **additive feature** $F: \mathcal{L} \to \mathbb{R}$ is defined additively based on a set of atomic contributions such as, for all $w \in \mathcal{L}$ generated by a unique parse tree $p_w$, one has

$$F(w) = \sum_{r \in p_w} f(r).$$

Features naturally include the number of occurrences $|w|_t$ of a given letter $t \in \Sigma$ in a word $w$. More complex features may require the design of an *ad hoc* grammar to accommodate highly-specific atomic contributions.

This allows to define the **bivariate generating function** of a language with respect to a
feature.

**Definition 5.1.6 (Bivariate generating function):** Given a language $\mathcal{L}$, the bivariate generating function $S(z, u)$ for a feature $F$ is

$$S(z, u) = \sum_{w \in \mathcal{L}} u^{F(w)} z^n = \sum_{n \ge 0} \sum_{k \ge 0} s_{n,k} u^k z^n \quad (5.7)$$

where $s_{n,k}$ is the number of words of length $n$ having value $k$ for the feature $F$.

The function $S(z, u)$ is solution of a system of equations similar to the constructs of Table 5.1, but with additional **monomials** $u^{f_r}$ being assigned to each production rule $r \in \mathcal{R}$ (unless $f_r = 0$, since $u^{f_r} = u^0 = 1$). Namely, for any non-terminal $S \in \mathcal{N}$, the system will feature an equation:

$$S(z, u) = \sum_{r: S \to w} u^{f_r} \prod_{x \in w} \begin{cases} T(z, u) & \text{if } x := T \in \mathcal{N}, \\ z & \text{if } x \in \Sigma, \\ 1 & \text{otherwise.} \end{cases}$$

Atomic contributions will accumulate at the exponent of $u$, ultimately leading to the term $u^{F(w)}$ in Equation (5.7).

**Expected value of a feature.** Solving the system, one obtains an expression for $S(z, u)$, from which information can be derived for the distribution of $F(W)$, the feature value of a random word. In particular, the **expected value** of $F$ for a uniformly distributed word in $\mathcal{L}_n$ is given by

$$\frac{[z^n] \frac{\partial S(z, u)}{\partial u} \Big|_{u \to 1}}{[z^n] S(z, 1)} = \frac{[z^n] \sum_{w \in \mathcal{L}} F(w) u^{F(w)-1} z^n \Big|_{u \to 1}}{[z^n] \sum_{w \in \mathcal{L}} z^n} = \frac{\sum_{w \in \mathcal{L}_n} F(w)}{s_n} \\ = \sum_{w \in \mathcal{L}_n} F(w) \times \mathbb{P}(W = w | n) = \mathbb{E}(F(W) | n)$$

Interestingly, algebraic functions remain algebraic upon partial derivatives, so singularity analysis can be used to obtain an asymptotic equivalent for the coefficients of the numerator.

### Example 22: Bivariate Motzkin words

To illustrate how bivariate generating functions can be used to analyze the content of random words, we consider the classic example of Motzkin words (see also Fig. 5.1), generated by the grammar

$$S \to (S) S | \bullet S | \varepsilon$$

Notice that, while ( and ) play a symmetrical role, the letter $\bullet$ seems quite different. Indeed, $\bullet$ is essentially responsible for the extension of sequences, while ( and ) emulate a parenthesis system, inducing tree-like structures. It is thus natural to pursue the *expected number of occurrences of* $\bullet$ *in a*
uniformly-distributed word of size n.

First, we define the atomic features as the number of occurrences of ● produced by each derivation:

$$r_1 : S \to (S) S \Rightarrow f_{r_1} = 0 \qquad r_2 : S \to \bullet S \Rightarrow f_{r_2} = 1 \qquad r_3 : S \to \varepsilon \Rightarrow f_{r_3} = 0$$

Then we translate the grammar into a system of functional equations involving by the bivariate generating function

$$S(z, u) = z^2 u^0 S(z, u)^2 + z u^1 S(z, u) + u^0 = z^2 S(z, u)^2 + z u S(z, u) + 1$$

Solving the system gives

$$S(z, u) = \sum_{n \ge 0} \sum_{k \ge 0} s_{n,k} u^k z^n = \frac{1 - zu - \sqrt{(1 - uz - 2z)(1 - uz + 2z)}}{2z^2}$$

One then derives by u, and set u = 1, to obtain the generating function counting the total number of occurrences of ● within Motzkin words of length n.

$$\begin{aligned} \frac{\partial S(z, u)}{\partial u} &= -\frac{z}{2z^2} + \frac{\sqrt{1 - uz - 2z}}{4z\sqrt{1 - uz + 2z}} + \frac{\sqrt{1 - uz + 2z}}{4z\sqrt{1 - uz - 2z}} \\ \left. \frac{\partial S(z, u)}{\partial u} \right|_{u \to 1} &= -\frac{z}{2z^2} + \frac{\sqrt{1 - 3z}}{4z\sqrt{1 + z}} + \frac{\sqrt{1 + z}}{4z\sqrt{1 - 3z}} \\ &= -\frac{z}{2z^2} + \frac{(1 - 3z) + (1 + z)}{4z\sqrt{1 + z} \times \sqrt{1 - 3z}} = -\frac{z}{2z^2} + \frac{1 - z}{2z\sqrt{1 + z}} \times \frac{1}{\sqrt{1 - 3z}} \end{aligned}$$

The resulting (mono variate) generating function features two singularities at $z = 1/3$ and $z = -1$ respectively, so the dominant one is $\rho = 1/3$. Applying Theorem 12 with $\alpha = 1/2$ and $U(z) = \sqrt{1 + z}/2z$ gives the following equivalent

$$[z^n] \left. \frac{\partial S(z, u)}{\partial u} \right|_{u \to 1} = \frac{1 - \rho}{2\rho\sqrt{1 + \rho\sqrt{\pi}}} \frac{\rho^{-n}}{\sqrt{n}} (1 + O(1/n))$$

On the other hand, for the denominator we have

$$S(z, 1) = S(z) = \frac{1 - z - \sqrt{(1 - 3z)(1 + z)}}{2z^2} \Rightarrow [z^n] S(z) = \frac{\sqrt{1 + \rho} \rho^{-n}}{4\rho^2 \sqrt{\pi} n \sqrt{n}} (1 + O(1/n))$$

Taking the ratio finally gives the expected number of occurrences of the letter ●

$$\begin{aligned} \mathbb{E}(|W| \cdot |n|) &= \frac{[z^n] \left. \frac{\partial S(z, u)}{\partial u} \right|_{u \to 1}}{[z^n] S(z)} = \frac{\frac{1-\rho}{2\rho\sqrt{1+\rho\sqrt{\pi}}} \frac{\rho^{-n}}{\sqrt{n}} (1+O(1/n))}{\frac{\sqrt{1+\rho}}{4\rho^2\sqrt{\pi}} \frac{\rho^{-n}}{n\sqrt{n}} (1+O(1/n))} \\ &= n \frac{2\rho(1-\rho)}{1+\rho} (1+O(1/n)) = \frac{n}{3} (1+O(1/n)) \end{aligned}$$

Quite surprisingly, the asymptotic proportion of ● is exactly 1/3, and so are those of ( and ) (having equal number of occurrences), despite the highly asymmetrical roles played by the three letters.

**Remark 5.1.3 (Higher Order Moments):** By iterating the partial derivative (compensating with u each time), **higher moments** can be extracted in a similar manner, granting access to the **variance**, **skewness**, **kurtosis** of the feature distribution induced by the uniform distribution over words of length n.
<figure><img src="image_42.png" /><figcaption>Figure 5.1: **Strong connectivity can be a matter of choice, as illustrated by two grammars for Motzkin words.** The top grammar has strongly-connected dependency graph, while the bottom grammar, although equivalent, induces a dependency graph which is not strongly-connected.</figcaption></figure>

### 5.1.4.2 Simple-type grammars

Shockingly, it is sometimes possible to skip steps 2, 3 and 4 of the symbolic method, and determine the type of singularity directly from properties of the grammar. To achieve *superlative laziness*, one first considers the dependencies in the grammar.

**Definition 5.1.7 (CFG Dependency Graph):** The **dependency graph** of a (weighted) grammar $G = (\mathcal{N}, \Sigma, \mathcal{R}, S_0)$ is a directed graph, whose vertices are non-terminals $(\mathcal{N})$, and whose arcs are any $S \to T$ such that $(S \to w) \in \mathcal{R}, T \in w$.

A critical phenomenon occurs when the dependency graph is **strongly-connected**, i.e. when any vertex can be reached from any other vertex through a sequence of arcs. Note that the absence of strong connectivity for a dependency graph reflects on the chosen grammar, and does not necessarily represent an intrinsic property of the language $\mathcal{L}$. It is thus the responsibility of the modeler to find a suitable grammar, as shown in example described in Figure 5.1.

Intuitively, the combinatorial explosion can be understood as **flowing** through the dependency graph, and the strong-connectivity implies the existence of a steady-state, where the dominant growth homogeneously contaminates all vertices (a.k.a. non-terminals). Under this condition, the growth of the language can be restricted to two possible categories, depending on whether the grammar is **linear**, i.e. at most one non-terminal in its production rules, or **tree-like**, i.e. at least one production rule with two or more produced non-terminal.
**Theorem 13 (Strongly-connected Asymptotics: Perron-Frobenius (69, pp 343) + Drmota-Lalley-Woods (56, 112, 202))**: Given a CFG G = ($\mathcal{N}$, $\Sigma$, $\mathcal{R}$, $S_0$) with strongly-connected dependency graph, only the following two types of dominant singularities may occur:

* If G is *linear*, then the dominant singularity $\rho$ is a *simple pole*, i.e. there exists $U(z)$ and $T(z)$, both analytic at $z = \rho$, with $U(\rho) \neq 0$ such that:

$$S(z) = \frac{U(z)}{1 - z/\rho} + T(z) \Rightarrow [z^n] S(z) = U(\rho) \times \rho^{-n}(1 + O(1/n));$$

* If G is *tree-like*, then the dominant singularity $\rho$ has *square-root type*, i.e. there exists $U(z)$ and $T(z)$ analytic at $z = \rho$, $k := \lim_{n \to +\infty} U(\rho) \neq 0$, such that:

$$S(z) = -U(z) \times \sqrt{1 - z/\rho} + T(z) \Rightarrow [z^n] S(z) = \frac{\kappa}{2\sqrt{\pi}} \times \frac{\rho^{-n}}{n\sqrt{n}}(1 + O(1/n)).$$

Such *stereotyped asymptotics* are immensely useful to study languages with large grammars, for which human patience and symbolic calculus can both reach their respective limits. In this case, computing the first coefficients of the generating function, e.g. using the recurrences of Figure 5.1, would allow to estimate $\kappa$ and $\rho$ to a sufficient precision in most applied contexts.

### 5.1.4.3 Weighted languages and asymptotics

Another extension, crucial in applied contexts, considers the introduction of weights, and was explored in a collaboration with Alain Denise and Michel Termier [50].

**Definition 5.1.8 (Weighted grammar):** A **Weighted Language** is generated by a context-free grammar G = ($\mathcal{N}$, $\Sigma$, $\mathcal{R}$, $S_0$), coupled with a **weight vector** $\pi$, associating a weight $\pi_t$ to each terminal $t \in \Sigma$.

The **weight** $\pi(w)$ **of a word** $w \in \Sigma^*$ is defined by

$$\pi(w) = \prod_{i=1}^{|w|} \pi_{w_i}.$$

The overall weight of a language $\mathcal{L}$ is then typically the quantity of interest, encapsulated in a **weighted generating function**

$$S_{\pi}(z) = \sum_{w \in \mathcal{L}} \pi(w) z^n = \sum_{n \ge 0} \sum_{w \in \mathcal{L}_n} \pi(w) z^n = \sum_{n \ge 0} s_n z^n$$
where $s_n$ now represents the accumulated weight of all words of length n in $\mathcal{L}$. Weighted generating functions associated with unambiguous grammars are the solutions of an algebraic system of functional equations. Those can be obtained using a slightly modified version of the constructs in Table 5.1, where **terminals rules** now account for their respective weight:

$$ (\text{Ordinary gen. fun.}) S(z) = z \rightarrow (\text{Weighted gen. fun.}) S_{\pi}(z) = \pi_t \times z. \quad (5.8) $$

While $S_{\pi}(z)$ can be of interest in itself, e.g. as its coefficients represent **partition functions**, one can also adopt an alternative probabilistic view and think of the weights as inducing a Boltzmann-Gibbs distribution over the class of words of size n.

**Definition 5.1.9 (Weighted distribution):** Let $G = (\mathcal{N}, \Sigma, \mathcal{R}, S_0, \pi)$ be an unambiguous CFG, generating a language $\mathcal{L}$, and n be a length, the **weighted distribution** is defined such that:

$$ P_{\pi}(w | n) = \frac{\prod_{i=1}^{|w|} \pi_{w_i}}{[z^n] S_{\pi}(z)} = \frac{\pi(w)}{s_n} $$

where $s_n = [z^n] S_{\pi} = \sum_{w \in \mathcal{L}_n} \pi(w)$.

The association of weights to terminal symbols (or, equivalently, to individual productions, as we showed [149]) grants more flexibility to the modeler in the design of a specification, particularly when weights are non-rational reals. Examples abound when real-valued weights are needed, e.g. when weights are optimized to achieve a targeted composition, or when classic Boltzmann-Gibbs distributions of the form $P(w) \propto e^{-E_w/\beta}$ are being investigated.

Moreover, their associated generating functions remain the solutions of algebraic systems of equations, so the entire *armada* of analytic combinatorics can be leveraged to analyze their asymptotic behavior.

**Proposition 14 (Asymptotics of Weighted Languages):** The asymptotic equivalents of Theorems 10, 11, 12, and 13 hold for weighted generating functions.

**Example 23: Expected number of secondary structures per RNA (212)**

The number of secondary structures of length n [194] only represents a crude upper bound for the number of secondary structures compatible with a given RNA sequence w, $|w| = n$. Indeed, canonical base-pairs can only form between pairs of nucleotides in

$$ B = \{(A, U), (U, A), (G, C), (C, G), (G, U), (U, G)\}, $$

forbidding certain secondary structures to form for certain sequences.

While the precise number of such structures vary substantially depending on the actual sequence,
it is possible to obtain a more precise estimate by considering the **expected number of secondary structures**. To that purpose, one may consider, and enumerate, **pairs** (*w*, *S*) of **compatible sequence and structure**, i.e. such that, for any base pairs in *S* the associated nucleotides in *w* belong to *B*. In other words, we need to find a grammar that jointly generates the sequence and structure, using an extended alphabet

$$ \Sigma^B = \{ (A, (C, (G, (U, )A, )C, )G, )U, \bullet A, \bullet C, \bullet G, \bullet U \} $$

leading to the grammar

$$ S \rightarrow (A S^{\ge 1})_U S \mid (U S^{\ge 1})_A S \mid (G S^{\ge 1})_C S \mid (C S^{\ge 1})_G S \\ | (G S^{\ge 1})_U S \mid (U S^{\ge 1})_G S \mid \bullet_A S \mid \bullet_C S \mid \bullet_G S \mid \bullet_U S \mid \varepsilon $$

Equivalently, one can simply consider the usual grammar of Example 21, and use weights $\pi_1(= 6, \pi_2) = 1, \pi_3 = 4$ to simulate the multiplicities of production rules. The associated system is then

$$ S_{\pi}(z) = \pi_{1}(z) z^{2} S_{\pi}(z) (S_{\pi}(z)^{\ge 1} - 1) + \pi_{2}\pi_{3}(z) + 1 = 6 z^{2} S_{\pi}(z) (S_{\pi}(z)^{\ge 1} - 1) + 4 S_{\pi}(z) + 1 $$

Solving the system gives the weighted generating function

$$ S_{\pi}(z) = \frac{1 - 4z + 6z^2 - \sqrt{1 - 8z + 4z^2 - 48z^3 + 36z^4}}{12z^2} $$

Singularity analysis can be performed, allowing estimates for the total number of compatible sequence/structure pairs

$$ S_{\pi}(z) = \frac{\kappa \rho^{-n}}{n\sqrt{n}} (1 + O(1/n)) $$

with $1/\rho \approx 8.164$. Dividing by the number of RNA sequences of size *n* gives the average number of structures compatible with a uniformly distributed sequence of length *n*, and we get

$$ E(\#\text{Compatible structures} | n) \in \Theta\left(\frac{2.04^n}{n\sqrt{n}}\right) $$

This approach can pushed even further, thanks to the flexibility of weighted/probabilistic generating functions, to produce asymptotics under more general random sequence models, including Bernoulli processes [212], or Markov chains [146].

## 5.2 Asymptotic combinatorics of RNA secondary structures

As mentioned in our introduction, RNA secondary structures have been the object of multiple works focusing on their asymptotic properties. We remind their formal definition, also defined in Section 1.2.3, for the sake of completeness.

**Definition 5.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**, $\forall(i, j) \neq (i', j') \in S: \{i, j\} \cap \{i', j'\} = \emptyset$;

* **Minimal distance** $\theta$ between paired nucleotides, $\forall(i, j) \in S: j - i > \theta$;

* **No pseudoknot** allowed, $\forall(i, j), (i', j') \in S, i < i': (j' < j)$ or $(j < i')$.

In the $\theta = 1$ case, i.e. no pairing of consecutive pairs, we obtain a popular combinatorial
<figure><img src="image_43.png" /><figcaption>Figure 5.2: **Decomposition of secondary structures and associated context-free grammar.** Note that one can always systematically rewrite a grammar to enforce cardinality constraints such as the one found in $S^{\ge 1}$.</figcaption></figure>

class (OEIS:A004148), admitting a wealth of equivalent characterizations, including:

1. Well-parenthesized expressions over $\{(\,,\,)\,,\bullet\}$ without occurrences of $\{\}$;

2. Positive walks from 0 to 0 taking steps $\{\uparrow, \downarrow, \to\}$ without peaks $(\uparrow\downarrow)$;

3. Unary-binary trees without **sherries** (simple binary tree, occurring as a leaf). . .

In particular, the representation as well parenthesized expressions immediately suggests a **decomposition for secondary structures**, ultimately leading to a grammar. Namely, consider a non-$\epsilon$ secondary structure, and focus on its first nucleotide b:

* If b is unpaired, then what follows is itself secondary structure;

* If b is paired to some partner b', then b and b' delimit an inner interval where a (non empty) structure can be found. Following the interval $[b, b']$ is another (unconstrained) secondary structure.

It follows that secondary structures can be split into two subsets, themselves built from (smaller) structures. Those subsets are disjoint since a given position cannot be simultaneously paired and unpaired, leading to an **unambiguous** decomposition. Their union, complemented with the empty structure $\epsilon$, is also **complete** with respect to secondary structures.

**Theorem 15 (Asymptotics of secondary structures (194))**: The number $s_n$ of secondary structures of length $n$ is asymptotically equivalent to

$$s_n = \sqrt{\frac{15 + 7\sqrt{5}}{8\pi}} \cdot \frac{\phi^n}{n\sqrt{n}} (1 + O(1/n)) \text{ where } \phi = \rho^{-1} \approx 2.62\dots$$

An illustrative proof of this result, using modern techniques of analytic combinatorics, is given in Example 21.
### 5.2.1 Expected 5'-3' distance (38)

A first application of the above principles considers the **5'-3' distance**, i.e. the graph distance induced by a(n ensemble of) secondary structure(s) between the ends of an RNA.

**Definition 5.2.2 (5'-3' distance):** Given an RNA secondary structure S, the **5'-3' distance**, graph distance induced by a(n ensemble of) secondary structure(s) between the ends of an RNA

Yoffe et al [207] argued that a small distance, independent of RNA length, between the two ends of RNA single-stranded molecules. The authors claimed that such a small distance was intrinsic to linear polymers, a claim supported by a heuristic argument arising from polymer physics, and empirically from the systematic MFE folding of randomly generated RNAs under various compositions. Such a short distance is indeed deemed crucial to the effective circularization of Eukaryotic messengers RNAs [196] through RNA-proteins and protein-protein complexes.

In collaboration with Peter Clote and Jean-Marc Steyaert [38], we decided to revisit this claim at three different levels: Theoretically, by checking its asymptotic validity in the homopolymer model; Algorithmically, by introducing a $\Theta(n^3)$ dynamic programming algorithm to compute the expected 5'-3' distance for a given sequence at the thermodynamic equilibrium; Empirically, by computing the graph distance on a collection of RNAs of documented structures.

Our asymptotic analysis relies on the design of a grammar which, by identifying nucleotides on the 5' to 3' shortest path, allowed the derivation of a bivariate generating function for the 5'-3' distance. We achieved this by a classic, if technical, application of the symbolic method. Essentially, we **duplicated the production rules**, distinguishing rules that generate the **exterior bases** from those whose content is nested within some base pair, and are as such not relevant to the expected 5'-3' distance. A similar reasoning allowed to capture general values for the minimum distance $\theta$ between two paired positions, and we obtained the following grammar:

$$S_{\theta} \rightarrow [R_{\theta}] S_{\theta} | \circ S_{\theta} | \varepsilon; \quad T_{\theta} \rightarrow (R_{\theta}) T_{\theta} | \bullet T_{\theta} | \varepsilon; \quad R_{\theta} \rightarrow (R_{\theta}) T_{\theta} | \bullet R_{\theta} | \bullet^{\theta} \quad (5.9)$$

Translating into a system of bivariate functional equation, one obtains

$$ \begin{aligned} S_{\theta}(z, u) &= z^2 u^2 R_{\theta}(z) S_{\theta}(z, u) + zuS_{\theta}(z, u) + 1 \\ T_{\theta}(z) &= z^2 R_{\theta}(z) T_{\theta}(z) + zT_{\theta}(z) + z^{\theta} \\ R_{\theta}(z) &= z^2 R_{\theta}(z) T_{\theta}(z) + zR_{\theta}(z) + z^{\theta}. \end{aligned} $$
Solving the system gives the bivariate generating function

$$S_{\theta}(z, u) = \sum_{n \ge 0} \sum_{k \ge 0} s_{k,n;\theta} u^k z^n = \frac{1}{1 - zu - \frac{u^2(1-2z+z^{\theta+2}-\sqrt{\Delta_{\theta}})}{2(1-z)}} \quad (5.10)$$

with $s_{k,n;\theta}$ the number of secondary structures of length n, with minimum base pair distance $\theta$, and having k nucleotides on the exterior face, and

$$\Delta_{\theta} = 1 - 4z + 4z^2 - 2z^{\theta+2} + 4z^{\theta+3} - 4z^{\theta+4} + z^{2\theta+4}. \quad (5.11)$$

In particular, one gets

$$S_{\theta}(z, 1) = \frac{1 - 2z + 2z^2 + 2z^{\theta} - z^{\theta+2} - \sqrt{\Delta_{\theta}}}{(1 - z)2z^2}.$$

Next, we use the classic bivariate analysis

$$E_{\theta}(z) = \left. \frac{\partial S_{\theta}(z, u)}{\partial u} \right|_{u=1} = \frac{P(z) - (2 - 5z + 4z^2 - 2z^{\theta+2} + z^{\theta+3})\sqrt{\Delta_{\theta}}}{2(1 - z)^2 z^4}$$

such that $e_n := [z^n] E_\theta(z)$ is the accumulated number of exterior bases across all RNAs of length n, and $P(z)$ is a polynomial of degree $\theta + 4$, whose precise value bears no consequence on the asymptotics. Applying Theorem 12, one gets

$$[z^n] S_{\theta}(z, 1) \sim -\frac{1}{(1-\rho)2\rho^2} \times [z^n] \sqrt{\Delta_{\theta}} \quad (5.12)$$

$$[z^n] E_{\theta}(z) \sim -\frac{(2 - 5\rho + 4\rho^2 - 2\rho^{\theta+2} + \rho^{\theta+3})}{2(1 - \rho)^2 \rho^4} \times [z^n] \sqrt{\Delta_{\theta}}. \quad (5.13)$$

where $\rho$ is the shared dominant singularity of $S_\theta(z)$ and $E_\theta$, i.e. the root of $\Delta_\theta$ having smallest modulus (real by Pringsheim). We obtain the following result by taking the ratio of estimates (5.13) and (5.12) to obtain the expected number of bases on the exterior base, to which we subtract 1 to obtain the graph distance.

**Theorem 16 (Expected 5'-3' distance - Homopolymer model (38))**: The **expected 5'-3' distance** $D_n$ over all RNA secondary structures of length n, assuming minimum distance $\theta$ between paired nucleotides, is given by

$$D_n = \frac{[z^n] E_\theta(z)}{[z^n] S_\theta(z, 1)} - 1 = \frac{2 - 5\rho + 3\rho^2 + \rho^3 - 2\rho^{\theta+2} + \rho^{\theta+3}}{(1 - \rho)\rho^2} (1 + O(1/n)) \quad (5.14)$$

where $\rho$ is the smallest modulus root of

$$\Delta_{\theta} := 1 - 4z + 4z^2 - 2z^{\theta+2} + 4z^{\theta+3} - 4z^{\theta+4} + z^{2\theta+4}. \quad (5.15)$$

As suspected by Yoffe et al [207], the expected value is constant, small and stable to different values of $\theta$. Of particular interest are the following numerical approximations,
for $\rho$ and $D_n$ when $\theta = 1$ and 3:

$$ \theta = 1 \rightarrow \rho \approx 0.381, D_n \approx 5.47 $$

$$ \theta = 3 \rightarrow \rho \approx 0.436, D_n \approx 4.15 $$

The expected distance can also be analyzed in a refined model introduced by Zuker and Sankoff [212] (see Example 23 for details), which captures base pairing rules that prevent certain structures to be adopted by a given sequence. In this model, a **stickiness parameter $\sigma$** represents the probability that a sequence, randomly generated by a Bernoulli process, forms a given pair ($\sigma := \sum_{\text{Valid BP}_s(b,b')} p_b \times p_{b'}$). Incorporating the stickiness in our above analysis allows to derive the **expected 5'-3' distance induced by a random RNA sequence**.

**Theorem 17 (Expected 5'-3' distance - Random Sequence model (38))**: The **expected 5' - 3' distance $D_n^\sigma$** in a random sequence of length $n$, generated by a Bernoulli process (stickiness $\sigma$), under a min. base pair distance $\theta$, is given by

$$ D_n^\sigma = \frac{\Phi_{\theta,\sigma}(\rho) - 4(\rho - 1)^3 \sigma^2 \rho^4}{4(\rho - 1)^3 \sigma^2 \rho^4} (1 + \mathcal{O}(1/n)) \quad (5.16) $$

where $\Phi_{\theta,\sigma}(z)$ is a (relatively) simple polynomial

$$ \Phi_{\theta,\sigma}(z) := 4 - 14z - 18z^2 + (2\sigma - 10)z^3 + (2 - 2\sigma)z^4 - 4\sigma z^{\theta+2} + 6\sigma z^{\theta+3} - 2\sigma z^{\theta+4}, $$

and $\rho$ is the root of smallest modulus of

$$ 1 - 4z + (6 - 2\sigma)z^2 + 4(\sigma - 1)z^3 + (\sigma - 1)^2 z^4 - 2\sigma z^{\theta+2} + 4\sigma z^{\theta+3} - 2\sigma(1 + \sigma)z^{\theta+4} + \sigma^2 z^{2\theta+4}. $$

Again, the introduction of the stickiness does not change the nature of the expected distance, which remain a constant of $n$.

**Our conclusions [38].** Interestingly, while both the asymptotics and equilibrium analysis of RNA sequences appeared to support the hypothesis of Yoffe *et al* [207], our empirical analyses of the STRAND database [8], consisting of structured non-coding RNAs, indicated in a positive observed correlation between the 5'-3' distance and RNA length. This observation suggests that structural RNAs may be under selective pressure to have larger 5'-3' distance than expected from pure polymer theory, warranting further analyses.

### 5.2.2 RNA network properties (180)

First introduced by Watts *et al* [195], **Small-world networks** are defined as satisfying the following, informally-stated, properties:
* the **shortest path distance** between any two nodes is “small”, e.g. six degrees of separation between any two persons;

* the average **clustering coefficient** is large, e.g. friends of a person tend to be friends of each other.

Small-world networks are ubiquitous in biology, sociology, and information technology; examples abound, including the neural network of C. elegans [195], the gene co-expression in S. cerevisiae [186], protein folding networks [22, 170]. . . For additional examples, see the excellent review of Albert and Barabási [4]. In particular, Wuchty [204] observed a small-world property for the low energy RNA secondary structure network of E. coli phe-tRNA, and hypothesized this property to be crucial for the structural kinetics induced by tRNA modifications.

In collaboration with Defne Surujon and Peter Clote [180], we investigated asymptotic properties of the RNA network, the ensemble of all RNA secondary structures connected by a **move set**.

**Definition 5.2.3 (Move set):** Given an RNA length n, a **move set**

$$MS = \{M_1, M_2, \dots\}$$

is a collection of **moves**, i.e. functions $M: S_n \to \mathcal{P}(S_n)$.

For any structure $S \in S_n$, the execution of a **move** returns structures $M(S) \subset S_n$, each resulting from the application of a (small) perturbation to S. Popular moves [72], in the context of a structure S, include:

* **Addition** of a base pair (i', j') such that $S \cup \{(i', j')\}$ is a valid structure;

* **Removal** of a base pair (i, j) $\in S$;

* **Shift** of (i, j) $\in S$ into (i', j'), i = i' or j = j', only if resulting in a valid structure;

We considered two **move sets**:

* $MS1 := \{\text{Addition, Removal}\}$;

* $MS2 := \{\text{Addition, Removal, Shift}\}$

Notice that, while individual moves may not be reversible, $MS1$ and $MS2$ are both symmetric, and allow to define a notion of **neighborhood** crucial to the definition of a network.
**Definition 5.2.4 (RNA Network):** Given a length n and a symmetric move set $\mathcal{MS}$, the RNA Network is the undirected graph $G_{\mathcal{MS}} = (V, E)$, where $V = S_n$ are secondary structures of size n, and E contains pairs of structures $\{S, S'\}$ such that $S' \in M(S)$ for some move $M \in \mathcal{MS}$ ($\Rightarrow S \in M'(S)$).

We particularly focused on the expected degree and clustering coefficient of large RNA networks, and rigorously proved that the RNA network is asymptotically not small-world. Following [40], a family $\{S_n, n = 1, 2, 3, ...\}$ of graphs must satisfy certain properties to qualify as small-world.

**Definition 5.2.5 (Small-world property (40)):** A family $\{G_n | n = 1, 2, 3, ...\}$ of graphs is small-world if and only if the following conditions hold:

A The largest distance between any two nodes in $G_n$ is in $O(\log |V|)$;

B The average degree of $G_n$ is in $O(\log |V|)$;

C The global clustering coefficient $\mathcal{C}_g(G_n)$ is bounded away from zero.

We remind that the **global clustering coefficient** $\mathcal{C}_g(G)$ of a graph $G = (V, E)$ is defined by Newman *et al* [141, Equation (77)] as:

$$ \mathcal{C}_g(G) = \frac{3 \times \#\text{triangles}}{\#\text{connected triples}} \qquad (5.17) $$

where a **triangle** is a set $\{x, y, z\} \in V$ of nodes such that $\{(x, y), (y, z), (z, v)\} \subset V$, and a **(connected) triple** is a set $\{x, y, z\}$ of nodes, such that $\{(x, y), (y, z)\} \subset V$.

Clearly, families of RNA networks induced by $\mathcal{MS}_1$ and $\mathcal{MS}_2$ both satisfy condition A, since any structure can be reached from any other one in $\Theta(n)$ operations (e.g. $\le n/2$ removals + $\le n/2$ additions), while $|V| = |S_n| \in \Theta(\rho^{-n}/n\sqrt{n})$.

**Expected degree.** The expected degree of a secondary structure, the average number of distinct structures in the neighborhood of a structure S, obviously depends on the move set. For $\mathcal{MS}_1$, starting from a structure S, any base pair in S can be removed, but only **co-accessible positions** (i.e. belonging to the same loop) can be safely added to S. For $\mathcal{MS}_2$, shift moves allow more flexibility, allowing the binding of positions involved in two adjacent loops. In both cases, we can already conclude that Condition B holds, since the degree of any S is in $O(n^3)$, i.e. bounded by a polynomial in n, while $|V| = |S_n|$ grows exponentially on n.

Still, a more precise analysis of the expected degree sheds some light on the techniques used to study the clustering coefficients, so we elaborate its main steps. Essentially, we
first observe that, in any structure $S \in S_n$, any base pair can be removed, so the overall number of removal moves within $S_n$ is exactly characterized by the accumulated number of base pairs in $S_n$.

Moreover, since additions and removals are reciprocal operations, there exists a bijection between pairs of structure/removal and structure/addition. The overall degree over $S_n$ is then $2 \times b_n$, where $b_n$ is the overall number of base pairs in $S_n$. Its associated generating function can be obtained by defining the number of base pairs as a feature, as shown in Section 5.1.4.1. Alternatively, we built an *ad hoc* grammar that generates all possible secondary structures with all possible *distinguished* base pair (here with $\theta = 1$):

$$S^* \rightarrow [S^{\geq \theta}] S | (S^{\bullet, \geq \theta}) S | (S^{\geq \theta}) S^* | \bullet S^* \\ S \rightarrow (S^{\geq \theta}) S | \bullet S | \varepsilon$$

Executing the various steps of the symbolic method, followed by singularity analysis, one obtains an asymptotic equivalent for $b_n$. Dividing $2 \cdot b_n$ by $|S_n|$ gives the expected degree $\delta_n^{MS1}$ under $\theta = 3$.

**Theorem 18 (Expected degree of the RNA network (180))**: The **expected degree** $\delta_n^{MS}$ of a random uniform node in the RNA network $G_n$ of a move set $MS$ is such that:

$$\delta_n^{MS1} = \kappa \cdot n(1 + \mathcal{O}(1/n)) \quad \text{with} \quad \kappa \approx 0.473;$$

$$\delta_n^{MS2} = \kappa \cdot n(1 + \mathcal{O}(1/n)) \quad \text{with} \quad \kappa \approx 1.56.$$

To compute the equivalent of $b_n$ for the shift move only, we built a grammar based on a similar idea as above, encoding the position and outcome of the move using an extended alphabet ($\star$ for the fixed end of the base pair, and $\langle \rangle$ for the shifting end). We obtained the following grammar, provably unambiguous and complete with respect to shifting moves:

$$\begin{align*}
\hat{S} &\rightarrow \hat{S}\bullet | (\hat{S}) | S(\hat{S}) | \hat{S}(R) | \hat{T} \\
\hat{T} &\rightarrow \star R\rangle\rangle | S \star R\rangle\rangle | \star R\rangle S\rangle | S \star R\rangle S\rangle | \langle\langle R\star | S\langle\langle R\star \\
&\quad | \langle S\langle R\star | S\langle S\langle R\star | \langle R \star R\rangle | S\langle R \star R\rangle \\
S &\rightarrow \bullet | S\bullet | (R) | S(R) \quad R \rightarrow \theta | R\bullet | (R) | S(R) \\
\theta &\rightarrow \bullet^\theta
\end{align*}$$

The resolution of the associated system, in itself, justifies investing in a symbolic computation software, but can nevertheless be done, leading to the asymptotics summarized in Theorem 19.
**Clustering coefficient.** Finally, we address the validity of Condition C, based on the value of the clustering coefficient, the last missing piece to assess the small world nature of homopolymer-based RNA networks. This requires the computation of the numbers of triangles and connected triples, two features which we capture in a similar strategy as above, and designing two unambiguous grammars whose generated languages are in bijection with the triangles and triples respectively. The associated grammars, which we spare the reader out of pure humanity, were tediously designed to exhaust (not only its designers but also) all pairs of valid consecutive moves in all contexts, this while capturing the proximity (triangle) or not (triple) of the initial and final structures.

**Theorem 19 (Triangles, triples and clustering in RNA networks (180))**: Under the $MS_2$ move set with $\theta = 3$, the RNA network $G_n$ is such that:

$$ \#triangles(G_n) = \kappa n (1 + O(1/n)) \quad \text{with } \kappa \approx 1.22 $$

$$ \#connected \ triples(G_n) = \lambda n^2 (1 + O(1/n)) \quad \text{with } \lambda \approx 0.112 $$

$$ \mathcal{C}_g(G_n) = \frac{3 \times \#triangles}{\#connected \ triples} \in \Theta\left(\frac{1}{n}\right) $$

It follows that $\mathcal{C}_g(G_n) \to 0$ as $n \to \infty$, meaning that Condition C is not satisfied. This result is robust to changes of $\theta$.

**Corollary 20 (RNA networks are not small world):**

The family of RNA networks is not small world.

**Conclusions.** Our result that, in an homopolymer model, RNA networks are not small world seem to contradict the conclusions of Wuchty [204] based on empirical observations. However, it must be noted that Wuchty [204] restricted its observations to the most stable conformations (low-free energy). Moreover, this observed discrepancy between a sequence-free (homopolymer) model and Wuchty's observations based on a limited subset of RNAs may be reconciled by the existence of a specific selective pressure for this family.

### 5.2.3 RNA shapes: Abstract representations of structures (120)

**RNA shapes**, introduced by Giegerich et al [82], represent a hierarchy of abstract representations for the secondary structures of RNA. Each secondary structure is assigned a shape at every level of abstractions, so that shapes at any level represent **equivalence classes** for the folding space. By providing a hierarchy coarse-grained descriptors, they
<figure><img src="image_44.png" /><figcaption>Figure 5.3: **Example of RNA shapes constructions from a secondary structure.** While $\pi'$ shapes (most detailed) merely contract unpaired regions and helices, $\pi$ shapes (coarsest) only retain the underlying tree backbone (gray tree).</figcaption></figure>

enable a variety of analysis, including a comparison-free clustering of secondary structures [27], e.g. generated by stochastic sampling [51] within the Boltzmann ensemble of low energy; or the alignment-free detection of recurrent conformations across homologous sequences [155]; or the indexation of RNA families for a faster attribution [102].

However, while shapes can be computed at various levels of abstractions (see Figure 5.3), and in linear time, the computational tractability, and statistical power, of a shape-based analysis is conditioned by a reasonable growth of the number of shapes populated by the method. To anticipate such a growth, in collaboration with Peter Cote and Andy Lorenz [120], we investigated the enumerative properties of RNA shapes at two extreme levels of abstractions, namely considering $\pi$ (coarsest) and $\pi'$ (most detailed).

**$\pi$-shapes.** As illustrated in Figure 5.3, $\pi$ shapes quite rudely disregard unpaired nucleotides, and contract consecutive helices separated only by an internal loop/bulge. They can then be seen as a well-parenthesized expression, forbidding directly nested (a.k.a. stacking) pairs of parentheses as these would be contracted by the abstraction process. RNA $\pi$ shapes (except for the shape of size 0, denoted as $\varepsilon$) can therefore be generated by the following grammar

$$S^{\pi} \rightarrow [T] S^{\pi} | [T] \quad T \rightarrow [T] S^{\pi} | \varepsilon$$

Solving the associated system, ultimately adding 1 to recover $\varepsilon$, gives

$$S^{\pi}(z) = \sum_{n \ge 0} s_{n}^{\pi} z^{n} = \frac{1 - z^{2} - \sqrt{1 - 2z^{2} - 3z^{4}}}{2z^{2}}$$

where $s_n^\pi$ is the total number of $\pi$-shapes over n characters, from which one draw two immediate conclusions.
<figure><img src="image_45.png" /><figcaption>Figure 5.4: Illustration of bijection between π shapes and Motzkin words. π shapes of size 2n + 2 can be interpreted as Dyck words of same length, and be turned into a binary tree, constrained so that the right child of an internal node cannot be a leaf. Removing leaves, one gets a partial binary tree with n edges, whose internal nodes must have a right child. Traversing the tree in preorder, emitting a letter ( when descending in a right child for the first time, a letter . for consecutive right descents, and a letter ) when going back up one level, one gets a Motzkin word (OEIS:A001006) of length n.</figcaption></figure>

**Theorem 21 (π-Shapes (120))**: The number $s_{2n}^{\pi}$ of π shapes over 2n characters obeys

$$s_{2n}^{\pi} = \frac{\sqrt{3}}{2\sqrt{\pi}} \cdot \frac{3^n}{n\sqrt{n}}(1 + \mathcal{O}(1/n))$$

Moreover, π-shapes of size 2n + 2 are in bijection with Motzkin words of length n.

The bijection can be made explicit, and we get:

$$\psi : [A] B \to \phi(A) \bullet \psi(B) \quad (B \neq \emptyset) \\ [A] \to \phi(A)$$

$$\phi : [A] B \to \phi(A) (\psi(B)) \\ \varepsilon \to \varepsilon.$$

The associated transformation can perhaps be better expressed as a traversal of a tree, as illustrated in by Figure 5.4.

**π'-shapes.** At the other end of the shape hierarchy, π' shapes are obtained from secondary structures by collapsing consecutive unpaired positions into a single dedicated character_. The language of π' shapes is then essentially a well-parenthesized expression with an neutral character, i.e. a Motzkin word (OEIS:A001006), excluding:

* Consecutive occurrences _ of unpaired regions;

* Nested/stacking pairs of parentheses [[ ... ]].

Expanding the classic Motzkin grammar to capture such constraints gives:

$$S \to U[T]S|U \quad T \to U[T]U[T]S|_{[T]} U_{[T]} U_{[T]} U_{[T]} \varepsilon \quad U \to \varepsilon$$
Solving the generating function, one gets

$$S^{\pi'}(z) = \frac{A(z) - B(z)\sqrt{\Delta}}{C(z)} \quad \text{with} \quad \Delta := 1 - 2z^2 - 6z^3 - 5z^4 - 6z^5 - z^6 + 2z^7 + z^8,$$

and A(z), B(z) and C(z) are polynomials, non-null at z = ρ the dominant singularity (also smallest root of Δ).

**Theorem 22 (π'-Shapes (120)):** The number $s_n^{\pi'}$ of π' shapes over n characters obeys

$$s_n^{\pi'} = \kappa \frac{\gamma^n}{n\sqrt{n}}(1 + O(1/n)) \text{ where } \gamma := \frac{1}{\rho} \approx 2.41 \text{ and } \kappa := \lim_{z \to \rho} \frac{B(z)}{2\sqrt{\pi C(z)}} \approx 0.985.$$

The chosen level of abstractions for shapes has therefore a strong impact on the number of representatives of the RNA folding landscapes. At the coarsest level, the number of π-shapes only grows in $Θ(1.73^n/n\sqrt{n})$, while the more informative π'-shapes induce a combinatorial explosion in $Θ(2.41^n/n\sqrt{n})$.

**Refined estimates for shapes compatible with a sequence.** The total number of shapes over n characters represents a very crude upper bound for the expected number of shapes returned by a predictive method when executed on an RNA consisting of n nucleotides. Indeed, certain sequences may not support the helices required to build a given shape. Ideally, one would therefore like to compute the **expected number of shapes** compatible with a sequence of length n. However, this turned out to be very challenging, due to apparently intrinsic difficulties to establish a suitable specification.

Indeed, computing an asymptotic equivalent for the expected number of shapes seems beyond the current reach of the symbolic method. For instance, contrary to RNA secondary structures, it does not seem possible to formulate an unambiguous grammars that enumerates all compatible shapes/sequences pairs. This is due to the fact that several secondary structures, all compatible with a given sequence and substantially different, may share a common shape at one or several levels, as shown below:

ACUACAGUGGUAGUACUUUAGAAUGUCUUAGA

$$S_1 = ((((..................)))) . ((..................)) . ((..................)) . \rightarrow \pi(S_1) = \text{[] []} \quad \pi(S_1) = \text{[] [] []}$$

$$S_2 = ((..................)) . ((..................)) . .................... \rightarrow \pi(S_2) = \text{[] []} \quad \pi(S_2) = \text{[] []}$$

$$S_3 = (((((..................)))) . (((((..................))))). \rightarrow \pi(S_3) = \text{[] []} \quad \pi(S_3) = \text{[] []}$$

More modestly, it is possible to refine our initial upper bound, by considering the **minimum length of a sequence supporting a given shape**. For instance, hairpin loops, identified by motifs [ ], are subject to a minimum distance θ, and thus can be thought to consume θ + 2 nucleotides. We can capture this phenomenon by replacing each occurrence
of $\epsilon$ in the above grammars, indicating the occurrence of a hairpin loop by $\bullet^\theta$. Summing over shapes associated with minimum size smaller than $n$, one obtains much refined estimates.

**Theorem 23 (Expected number of shapes (120))**: The expected number $e_n^\pi$ (resp. $e_n^{\pi'}$) of $\pi$ shapes (resp. $\pi'$ shapes), compatible with an RNA of length $n$ under a min base pairing distance $\theta = 3$, obeys:

$$e_n^\pi \le \frac{1.28 \times 1.81^n}{n\sqrt{n}}(1 + \mathcal{O}(1/n))$$

$$e_n^{\pi'} \le \frac{2.44 \times 1.32^n}{n\sqrt{n}}(1 + \mathcal{O}(1/n))$$

While these upper bound are considerably tighter than the total number of shapes of length $n$, they still appear overly pessimistic in comparison to the empirical observations of Voss et al [190]. In this work, the authors conjectured growths in $\Theta(1.1^n)$ and $\Theta(1.16^n)$ for $\pi$ and $\pi'$ shapes respectively, motivating the development of new analysis techniques for this pragmatic class of conformational descriptors.

### 5.2.4 Enumerating designable structures (206)

The **inverse folding** of RNA, i.e. the computational search for an RNA sequence that folds stably and preferentially into a **target structure**, is a classic test piece in RNA bioinformatics. Its motivation arises from applications ranging from synthetic biology [182] to RNA therapeutics [203] through systems biology [65] and nanotechnologies [83]. It is NP-hard even under very simple energy models [20], and is typically approached at a practical level through a large collection of methods, based on metaheuristics, constraint programming or random generation paradigms [34].

Since the overarching goal of RNA design is to produce functional molecules, it is tempting to adopt an inverted perspective on the RNA design, and treat the limitations induced by popular sequence generation models as constraints weighing on evolution itself. Indeed, the sequence/structure relationship in RNA is used as a metaphor for the genotype/phenotype map in theoretical evolutionary studies [105] (neutral network theory). In this context, the total [194] or expected [212] number of secondary structures is used as a baseline to quantify the **phenotype space**, the set of functions achieved by a gene/molecule. However, this set is largely overestimated since many structures cannot be adopted as a functional structure, e.g. due to their intrinsic unstability.

In collaboration with Hua-Ting Yao, Cédric Chauve and Mireille Régnier [206], we addressed the quantification of **designable structures**, i.e. structures for which no RNA sequence is deemed acceptable which respect to the objectives of inverse folding. More specifically, following the seminal work of Dirks et al [54], notions of **defects** are used as
objectives for design. Defect metrics quantify, for a given candidate sequence, the existence of competitive alternative structures to the target structure. Examples include the **free-energy defect**, the free-energy difference, for a candidate sequence, between the target structure and its best competitor; or the **probably defect**, the overall Bolzmann probability of not adopting the target structure at the thermodynamic equilibrium. This allows to define a **design criterion**, which only accepts sequences whose defect does not exceed some predefined **tolerance** $\epsilon > 0$ for the defect $\mathcal{D}$.

**Definition 5.2.6 (Designable secondary structure):** Given a defect $\mathcal{D}$ and a tolerance $\epsilon > 0$, a secondary structure $S^* \in S_n$ is **designable** if and only if there exists $w \in \Sigma^n$ such that $\mathcal{D}(w, S^*) < \epsilon$.

Since their computation requires the joint consideration of: i) a realistic energy model [185]; and ii) the entire space of alternative structures, the computation of defects usually involves the execution of a dynamic programming algorithm. Due to their reliance on dynamic programming, all modern notions of defect are **monotonous with respect to loops**. Namely, as soon as a sequence/structure pair exceeds a certain **defect tolerance** within a (combination of) loop(s), coined a **local obstruction**, it cannot be extended into a larger structure sequence that would meet tolerance. This fact was already observed by Aguirre-Hernández et al [1], identifying two local obstructions for the free-energy defect, and by Hales et al [86] in a base pair maximization model. We verified the **ubiquity of obstructions**, proposing a *brute force* algorithm running in $\Theta(\Phi^k 4^k k^3)$, $\Phi = (1 + \sqrt{5})/2$ (a.k.a. the golden ratio), for computing the intrinsic defect of all (combinations of) loops over k nucleotides, ultimately resulting in a list $\mathcal{F}$ of local obstructions for a given defect and threshold.

**Remark 5.2.1:** Regardless of the chosen defect, the existence of a polynomial time algorithm for the enumeration all local obstructions would be highly surprising. Indeed, the existence of a complete list of obstructions of size n, having cardinality polynomial on n, would suggest a polynomial-time algorithm, using basic motif search in trees, for the decision version of the NP-hard inverse folding problem [20]. Possible alternatives, such as fixed-parameter algorithms, remain a possibility that could be worthy of being explored.

From a collection $\mathcal{F}$ of local obstructions, the next step is to consider the restricting effect of forbidden motifs on the structure/phenotype space. To that purpose, we adopted a tree representation, where the occurrence of a motif $m \in \mathcal{F}$ within a structure can be expressed as the exact match of a tree motif associated with m, as illustrated by Figure 5.5. It is well-known that motifs in trees asymptotically occur a number of time linear on the tree size, almost irrespectively of their precise definition, and that forbidding some motif induces a decay of the exponential growth [35, 39]. The following result ensues immediately.
<figure><img src="image_46.png" /><figcaption>Figure 5.5: Graph and tree representations of an RNA secondary structure S. Our notion of occurrence of a motif in a structure requires an exact match, such that the induced subgraph m₂ of a motif m₁ ∈ S does not automatically occur in S.</figcaption></figure>

**Theorem 24 (Exponentially low density of designable structures (206)):** For any defect D and tolerance ε, the proportion of **designable secondary structure** in Sⁿ is exponentially decreasing with n, i.e. there exists β < 1 such that

$$ \frac{|\{S \in S_n | S \text{ is designable}\}|}{|S_n|} < \beta^n. $$

To establish this result, the set of secondary structures that avoids any motif in F can be used, as it represents a superset for the set of designable structures. It is generated by the grammar:

$$ \begin{aligned} S &\rightarrow (T) S | \bullet S | \varepsilon \\ T &\rightarrow S \setminus F \end{aligned} $$

with $F := \{m' | \forall m \in F, m = (m')\}$. This translates into the following system:

$$ \begin{aligned} S(z) &= z^2 T(z) S(z) + z S(z) + 1 \\ T(z) &= S(z) - F(z, T) \end{aligned} $$

where $F(z, T)$ is the gen. fun. of enclosed structures within motifs in F, such that

$$ F(z, T) = \sum_{m \in \mathcal{F}} z^{\gamma(m)} T^{\delta(m)} - c(z, T) $$

where $\gamma(m)$ (resp. $\delta(m)$) is the size (resp. number of paired leaves) of the motif m', and $c(z, T)$ is a correcting term to account for potential overlaps. Solving the equation numerically gives the asymptotic equivalents listed in Table 5.2.
<table><caption>Table 5.2: Consequence of obstructions of size up to 12 on proportion of designable secondary structures. (*) Proportions of designable sequences computed using an assumption of equal constants for the asymptotic leading terms of the number of secondary structures, respectively allowing and forbidding local obstructions.</caption><thead><tr><th rowspan="2">Defect</th><th rowspan="2">ε</th><th rowspan="2">obstructions</th><th rowspan="2">ρ</th><th rowspan="2">Asymptotic<br>equivalent</th><th colspan="7">Proportion of designable structures* (upper bound)</th></tr><tr><th>Equiv.</th><th>P<sub>10</sub> (%)</th><th>P<sub>50</sub> (%)</th><th>P<sub>100</sub> (%)</th><th>P<sub>200</sub> (%)</th><th>P<sub>500</sub> (%)</th><th>P<sub>1000</sub> (%)</th></tr></thead><tbody><tr><td>Free-energy</td><td>1</td><td>104</td><td>0.449</td><td>Θ(<sup>2.226<sup>n</sup></sup>/<sub>n√n</sub>)</td><td>0.973<sup>n</sup></td><td>76.1</td><td>25.4</td><td>6.48</td><td>4.19·10<sup>-1</sup></td><td>1.14·10<sup>-4</sup></td><td>1.30·10<sup>-10</sup></td></tr><tr><td>Probability</td><td>.5</td><td>117</td><td>0.450</td><td>Θ(<sup>2.224<sup>n</sup></sup>/<sub>n√n</sub>)</td><td>0.972<sup>n</sup></td><td>75.3</td><td>24.2</td><td>5.84</td><td>3.41·10<sup>-1</sup></td><td>6.81·10<sup>-5</sup></td><td>4.64·10<sup>-11</sup></td></tr><tr><td>Probability</td><td>.1</td><td>152</td><td>0.460</td><td>Θ(<sup>2.176<sup>n</sup></sup>/<sub>n√n</sub>)</td><td>0.95<sup>n</sup></td><td>59.9</td><td>7.69</td><td>0.59</td><td>3.51·10<sup>-3</sup></td><td>7.27·10<sup>-10</sup></td><td>5.29·10<sup>-21</sup></td></tr><tr><td>Probability</td><td>.01</td><td>174</td><td>0.481</td><td>Θ(<sup>2.078<sup>n</sup></sup>/<sub>n√n</sub>)</td><td>0.908<sup>n</sup></td><td>38.1</td><td>0.80</td><td>6.44·10<sup>-3</sup></td><td>4.14·10<sup>-7</sup></td><td>1.10·10<sup>-19</sup></td><td>1.22·10<sup>-40</sup></td></tr></tbody></table>
# Chapter 6

## Conclusion and perspectives

In this Habilitation, I have attempted to provide a unified perspective over a series of contributions in ensemble dynamic programming, random generation, analysis and design of algorithms with applications in Bioinformatics, focusing on RNA folding from thermodynamic principles and RNA folding.

These research projects, spanning more than a decade of research, more than a dozen of distinct collaborations and 50+ manuscripts published in journals and selective conference, have allowed me to draw fruitful connections between my initial background in computer science and discrete applied mathematics, and my growing expertise and visibility in RNA Bioinformatics. Beyond their absolute value as scientific contributions and tools for the community, they have contributed to build a level of clarity in my objectives, and positioning with respect to my communities, allowing me to (more) confidently envision the future mentoring of young researchers.

**Current projects.** Over the next couple of years, my research will foreseeably be dedicated to two main projects in RNA bioinformatics, namely the **connection between design and RNA evolution** and the **integrative structure prediction informed by probing data**, under respective grants from the French Agence Nationale de la Recherche (DECRYPTED and PARNAssus projects) awarded for the two projects for the 2020–2024 period.

In the former, we will explore with Simona Cocco, Rémi Monasson and Bruno Sargueil to which extent modern methods in evolutionary analysis, based on the concept of **direct information**, can be used to recover and unravel design principles resulting from negative selective pressures. I also wish to revisit negative design as a null model for assessing the significance of observations, as mentioned in Chapter 3. This will require revisiting the controlled random generation of a NP-hard problem, hopefully helped by my recent enthusiasm for parameterized complexity algorithms.
The latter project, implemented in collaboration with Ronny Lorenz and Bruno Sargueil, attempts to define new methods to exploit information produced in **unconventional settings for RNA probing**. Those include probing experiments in different ionic conditions where pseudoknots are suspected (or new tests to assess their presence), produced using diverse reagents revealing different aspects of RNA architectures or acting with different speeds. One of the requirements for such an endeavor is the design of polynomial DP-based classes of pseudoknots, using conceptual tools introduced in Chapter 2, efficient and focusing on geometrically-feasible structures, enabling reasonable ensemble analyses.

At a more methodological level, I am increasingly convinced that an efficient pursue of
finalized projects in RNA Bioinformatics requires to limit the distraction represented by
the implementation/debugging of complex DP algorithms (by students having limited
experience and patience). For this reason, I plan to invest some of my time to the
development of a **generic DP framework** to allow the production of efficient low-level
code from a generic DP description. Existing frameworks are indeed either insufficiently
generic, or not easily amenable to support desirable features like random generation.
Such a development would also provide motivation to explore connections with formal
calculus, from which some of my contributions are already directly originating (through
enumerative combinatorics).

Finally, the recent development of methods combining sampling and clustering, follow-
ing Ding et al. [52] have led us to identify strong limitations in performances, experienced
when using both agglomerative and divisive algorithms for unsupervised machine. I
am therefore currently venturing in the *darkest pits* of machine learning, hoping to find
techniques (e.g. embeddings) that will allow an efficient processing of large sampled
sets of structure, hoping to increase the robustness, and reproducibility of downstream
analyses.

**Team development.** This decade since my PhD has finally allowed to grow as a scientist, from the postdoc-like status associated with a junior researcher CNRS position, to securing my own research fundings and supervising my student, and finally to taking the acting (2016) and permanent (2017–2018) scientific lead of the AMIBio team at LIX, culminating in the recent recruitment of Sebastian Will.
Bibliography

[1] R. Aguirre-Hernández, H. H. Hoos, and A. Condon. Computational RNA secondary structure design: empirical complexity and improved methods. BMC Bioinformatics, 8:34, 2007.

[2] A. V. Aho and M. J. Corasick. Efficient string matching: an aid to bibliographic search. Communications of the ACM, 18(6):333–340, jun 1975.

[3] T. Akutsu. Dynamic programming algorithms for RNA secondary structure prediction with pseudo-knots. Discrete Applied Mathematics, 104(1):45–62, 2000.

[4] R. Albert and A.-L. Barabási. Statistical mechanics of complex networks. Reviews of modern Physics, 74:47–97, 2002.

[5] D. Aldous. On the markov chain simulation method for uniform combinatorial distributions and simulated annealing. Probability in the Engineering and Informational Sciences, 1(1):33–46, 1987.

[6] S. F. Altschul and B. W. Erickson. Significance of nucleotide sequence alignments: a method for random sequence permutation that preserves dinucleotide and codon usage. Molecular biology and evolution, 2:526–538, Nov. 1985.

[7] M. Andronescu, A. P. Fejes, F. Hutter, H. H. Hoos, and A. Condon. A new algorithm for RNA secondary structure design. Journal of Molecular Biology, 336(3):607–624, 2004.

[8] M. Andronescu, V. Bereg, H. H. Hoos, and A. Condon. RNA strand: the RNA secondary structure and statistical analysis database. BMC bioinformatics, 9:340, Aug. 2008.

[9] A. Avihoo, A. Churkin, and D. Barash. RNAexinv: An extended inverse RNA folding from shape and physical attributes to sequences. BMC Bioinformatics, 12(1):319, 2011.

[10] J. K. Baker. Trainable grammars for speech recognition. The Journal of the Acoustical Society of America, 65(S1):S132–S132, 1979.

[11] C. Banderier and M. Drmota. Formulae and asymptotics for coefficients of algebraic functions. Combinatorics, Probability and Computing, 24(1):1–53, 2015.

[12] 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.

[13] E. Barcucci, R. Pinzani, and R. Sprugnoli. The random generation of directed animals. Theoretical Computer Science, 127(2):333–350, 1994.

[14] R. Bellman. The theory of dynamic programming. Bulletin of the American Mathematical Society, 60(6):503–515, 1954.

[15] M. Bendkowski, O. Bodini, and S. Dovgal. Polynomial tuning of multiparametric combinatorial samplers. In 2018 Proceedings of the Fifteenth Workshop on Analytic Algorithmics and Combinatorics (ANALCO), pages 92–106. SIAM, 2018.

[16] P. Berenbrink and T. Sauerwald. The weighted coupon collector's problem and applications. In 15th International Computing and Combinatorics Conference (COCOON'10), 2009.

[17] G. Blin, A. Denise, S. Dulucq, C. Herrbach, and H. Touzet. Alignments of RNA structures. IEEE/ACM Trans. Comput. Biology Bioinform., 7(2):309–322, 2010.

[18] 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.
[19] M. Bon and H. Orland. TT2NE: a novel algorithm to predict RNA secondary structures with pseudo-knots. Nucleic Acids Research, 39(14):e93-e93, 05 2011.

[20] É. Bonnet, P. Rzążewski, and F. Sikora. Designing RNA secondary structures is hard. In Research in Computational Molecular Biology 22nd Annual International Conference, RECOMB 2018, volume 10812 of Lecture Notes in Computer Science, pages 248-250, Paris, 2018. Springer.

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

[22] G. R. Bowman and V. S. Pande. Protein folded states are kinetic hubs. Proc. Natl. Acad. Sci. U.S.A., 107(24):10890-10895, June 2010.

[23] R. Bundschuh and R. Bruinsma. Melting of branched RNA molecules. Physical review letters, 100 (14):148101, 2008.

[24] R. Bundschuh and T. Hwa. RNA secondary structure formation: A solvable model of heteropolymer folding. Phys. Rev. Lett., 83:1479-1482, Aug 1999.

[25] R. Bundschuh and T. Hwa. Statistical mechanics of secondary structures formed by random RNA sequences. Physical review. E, Statistical, nonlinear, and soft matter physics, 65:031903, Mar. 2002.

[26] R. Bundschuh and T. Hwa. Phases of the secondary structures of RNA sequences. EPL (Europhysics Letters), 59(6):903, 2002.

[27] H. A. Burbano and E. Andrade. Analysis of tRNA abstract shapes of precursor/derivative amino acids in archaea. Gene, 396:75-83, July 2007.

[28] A. Busch and R. Backofen. INFO-RNA—a fast approach to inverse RNA folding. Bioinformatics, 22 (15):1823-31, 2006.

[29] S. Cao and S. Chen. Predicting RNA pseudoknot folding thermodynamics. Nucleic Acids Research, 34(9):2634-2652, 2006.

[30] 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.

[31] 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.

[32] 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.

[33] 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.

[34] 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.

[35] F. Chyzak, M. Drmota, T. Klausner, and G. Kok. The distribution of patterns in random trees. Combinatorics, Probability and Computing, 17(1):21-59, 2008.

[36] P. Clote. Combinatorics of saturated secondary structures of RNA. Journal of computational biology : a journal of computational molecular cell biology, 13:1640-1657, Nov. 2006.

[37] P. Clote, F. Ferré, E. Kranakis, and D. Krizanc. Structural RNA has lower folding energy than random RNA of the same dinucleotide frequency. RNA (New York, N.Y.), 11:578-591, May 2005.

[38] P. Clote, Y. Ponty, and J.-M. Steyaert. Expected distance between terminal nucleotides of RNA secondary structures. Journal of Mathematical Biology, 65(3):581-99, Sept. 2012.

[39] G. Collet, J. David, and A. Jacquot. Random sampling of ordered trees according to the number of occurrences of a pattern. Submitted, 2018.

[40] R. Cont and E. Tanimura. Small-world graphs: characterization and alternative constructions. Adv. in Appl. Probab., 40(4):939-965, 2008.

[41] J. Cooley and J. Tukey. An algorithm for the machine calculation of complex fourier series. Mathematics of Computation, 19(90):297-301, 1965.

[42] J. Cooley, P. Lewis, and P. Welch. The finite Fourier transform. IEEE Transactions on Audio and Electroacoustics, 17(2):77-85, June 1969.
[43] J. Cupal, I. L. Hofacker, and P. F. Stadler. Dynamic programming algorithm for the density of states of RNA secondary structures. In *Proceedings of the German Conference on Bioinformatics*, GCB, pages 184–186, Leipzig, Germany, Sept. 1996.

[44] K. Darty, A. Denise, and Y. Ponty. VARNA: Interactive drawing and editing of the RNA secondary structure. *Bioinformatics*, 25(15):1974–5, Aug. 2009.

[45] N. S. de Groot, A. Armaos, R. Graña-Montes, M. Alriquet, G. Calloni, R. M. Vabulas, and G. G. Tartaglia. RNA structure drives interaction with proteins. *Nature Communications*, 10(1), Jul 2019.

[46] 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.

[47] K. E. Deigan, T. W. Li, D. H. Mathews, and K. M. Weeks. Accurate SHAPE-directed RNA structure determination. *Proc Natl Acad Sci U S A*, 106(1):97–102, 2009.

[48] A. Denise. Génération aléatoire et uniforme de mots. *Discrete Mathematics*, 156:69–84, 1996.

[49] A. Denise and P. Zimmermann. Uniform random generation of decomposable structures using floating-point arithmetic. *Theoretical Computer Science*, 218(2):233–248, 1999.

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

[51] Y. Ding and C. E. Lawrence. A statistical sampling algorithm for RNA secondary structure prediction. *Nucleic acids research*, 31:7280–7301, Dec. 2003.

[52] Y. Ding, C. Y. Chan, and C. E. Lawrence. RNA secondary structure prediction by centroids in a boltzmann weighted ensemble. *RNA (New York, N.Y.)*, 11:1157–1166, Aug. 2005.

[53] R. Dirks and N. Pierce. A partition function algorithm for nucleic acid secondary structure including pseudoknots. *J Comput Chem*, 24:1664–1677, 2003.

[54] R. M. Dirks, M. Lin, E. Winfree, and N. A. Pierce. Paradigms for computational nucleic acid design. *Nucleic Acids Research*, 32(4):1392–1403, 2004.

[55] C. B. Do, D. A. Woods, and S. Batzoglou. Contrafold: RNA secondary structure prediction without physics-based models. *Bioinformatics (Oxford, England)*, 22:e90–e98, July 2006.

[56] M. Drmota. Systems of functional equations. *Random Structures and Algorithms*, 10(1-2):103–124, 1997.

[57] 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.

[58] 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.

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

[60] P. Duchon, P. Flajolet, G. Louchard, and G. Schaeffer. Boltzmann samplers for the random generation of combinatorial structures. *Combinatorics, Probability and Computing*, 13(4-5):577–625, 2004.

[61] J. Duraj. Random walks in cones: The case of nonzero drift. *Stochastic Processes and their Applications*, 124(4):1503–1518, apr 2014.

[62] S. R. Eddy and R. Durbin. RNA sequence analysis using covariance models. *Nucleic Acids Research*, 22(11):2079–2088, 1994.

[63] A. Esmaili-Taheri and M. Ganjtabesh. Erd: a fast and reliable tool for RNA design including constraints. *BMC bioinformatics*, 16:20, Jan. 2015.

[64] A. Esmaili-Taheri, M. Ganjtabesh, and M. Mohammad-Noori. Evolutionary solution for the RNA design problem. *Bioinformatics*, 30(9):1250–1258, 2014.

[65] S. Findeiß, M. Etzel, S. Will, M. Mörl, and P. F. Stadler. Design of artificial riboswitches as biosensors. *Sensors (Basel, Switzerland)*, 17(9):E1990, Aug. 2017.
[66] A. V. Finkelstein and M. A. Roytberg. Computation of biopolymers: a general approach to different problems. **Biosystems**, 30(1-3):1-19, 1993.

[67] P. Flajolet. Analytic models and ambiguity of context-free languages. **Theoretical Computer Science**, 49:283-309, 1987.

[68] P. Flajolet and A. M. Odlyzko. Singularity analysis of generating functions. **SIAM J. Discrete Math.**, 3(2):216-240, 1990.

[69] P. Flajolet and R. Sedgewick. Analytic combinatorics. cambridge University press, New York, NY, USA, 1 edition, 2009. ISBN 0521898064, 9780521898065.

[70] P. Flajolet, D. Gardy, and L. Thimonier. Birthday paradox, coupon collectors, caching algorithms and self-organizing search. **Discrete Applied Mathematics**, 39(3):207-229, 1992.

[71] P. Flajolet, P. Zimmermann, and B. V. Cutsem. A calculus for the random generation of labelled combinatorial structures. **Theoretical Computer Science**, 132(1):1-35, 1994.

[72] C. Flamm, W. Fontana, I. Hofacker, and P. Schuster. RNA folding at elementary step resolution. **RNA**, 6:325-338, 2000.

[73] C. Flamm, I. L. Hofacker, S. Maurer-Stroh, P. F. Stadler, and M. Zehl. Design of multistable RNA molecules. **RNA (New York, N.Y.)**, 7:254-265, Feb 2001.

[74] W. Fontana, P. F. Stadler, E. G. Bornberg-Bauer, T. Griesmacher, I. L. Hofacker, M. Tacker, P. Tarazona, E. D. Weinberger, and P. Schuster. RNA folding and combinatorial landscapes. **Phys. Rev. E**, 47: 2083-2099, Mar 1993.

[75] E. Freyhult, V. Moulton, and P. Clote. Boltzmann probability of RNA structural neighbors and riboswitch detection. **Bioinformatics (Oxford, England)**, 23:2054-2062, Aug. 2007.

[76] D. Gamerman and H. F. Lopes. Markov chain Monte Carlo: stochastic simulation for Bayesian inference. Chapman and Hall/CRC, 2006.

[77] J. A. Garcia-Martin, P. Clote, and I. Dotu. RNAiFOLD: a constraint programming algorithm for RNA inverse folding and molecular design. **Journal of Bioinformatics and Computational Biology**, 11(2): 1350001, 2013.

[78] 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.

[79] Q. Ge and D. Štefankovič. A graph polynomial for independent sets of bipartite graphs. **Combinatorics, Probability and Computing**, 21(05):695-714, 2012.

[80] R. Giegerich and C. Meyer. Algebraic dynamic programming. In International Conference on Algebraic Methodology and Software Technology, pages 349-364. Springer, 2002.

[81] R. Giegerich and H. Touzet. Modeling dynamic programming problems over sequences and trees with inverse coupled rewrite systems. **Algorithms**, 7(1):62-144, Mar 2014.

[82] R. Giegerich, B. Voss, and M. Rehmsmeier. Abstract shapes of RNA. **Nucleic acids research**, 32: 4843-4851, 2004.

[83] S. Grabbe, H. Haas, M. Diken, L. M. Kranz, P. Langguth, and U. Sahin. Translating nanoparticulate-personalized cancer vaccines into clinical applications: case study with RNA-lipoplexes for the treatment of melanoma. **Nanomedicine (London, England)**, 11:2723-2734, Oct. 2016.

[84] S. Griffiths-Jones, A. Bateman, M. Marshall, A. Khanna, and S. R. Eddy. RFAM: an RNA family database. **Nucleic Acids Research**, 31(1):439-441, 2003.

[85] 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.

[86] 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.

[87] M. Hamada, H. Kiryu, K. Sato, T. Mituyama, and K. Asai. Prediction of RNA secondary structure using generalized centroid estimators. **Bioinformatics (Oxford, England)**, 25:465-473, Feb. 2009.

[88] 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.
[89] S. Hammer, C. Günzel, M. Mörl, and S. Findeiß. Evolving methods for rational design of functional RNA molecules. **Methods**, 161:54–63, 2019.

[90] 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.

[91] M. Hellmuth, D. Merkle, and M. Middendorf. Extended shapes for the combinatorial design of RNA sequences. **International journal of computational biology and drug design**, 2:371–384, 2009.

[92] C. Herrbach, A. Denise, and S. Dulucq. Average complexity of the Jiang-Wang-Zhang pairwise tree alignment algorithm and of a RNA secondary structure alignment algorithm. **Theor. Comput. Sci.**, 411(26-28):2423–2432, 2010.

[93] P. G. Higgs and N. Lehman. The RNA world: molecular cooperation at the origins of life. **Nature Reviews Genetics**, 16(1):7–17, nov 2014.

[94] F. Hildebrand, A. Meyer, and A. Eyre-Walker. Evidence of selection upon genomic gc-content in bacteria. **PLoS genetics**, 6:e1001107, Sept. 2010.

[95] I. L. Hofacker, W. Fontana, P. F. Stadler, L. S. Bonhoeffer, M. Tacker, and P. Schuster. Fast folding and comparison of RNA secondary structures. **Monatshefte für Chemie/Chemical Monthly**, 125(2): 167–188, 1994.

[96] I. L. Hofacker, P. Schuster, and P. F. Stadler. Combinatorics of RNA secondary structures. **Discrete Applied Mathematics**, 88(1-3):207–237, 1998.

[97] C. Höner zu Siederdissen. Sneaking around concatmap: Efficient combinators for dynamic programming. **SIGPLAN Not.**, 47(9):215–226, Sept. 2012.

[98] L. Huang and D. Chiang. Better k-best parsing. In *Proceedings of the Ninth International Workshop on Parsing Technology*, pages 53–64. Association for Computational Linguistics, 2005.

[99] L. D. Hurst and A. R. Merchant. High guanine-cytosine content is not an adaptation to high temperature: a comparative analysis amongst prokaryotes. **Proceedings. Biological sciences**, 268:493–497, Mar. 2001.

[100] S. Ieong, M. yang Kao, T. wah Lam, W. kin Sung, and S. ming Yiu. Predicting RNA secondary structures with arbitrary pseudoknots by maximizing the number of stacking pairs. **Journal Of Computational Biology**, 10(6):981–995, 2003.

[101] 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.

[102] S. Janssen, J. Reeder, and R. Giegerich. Shape based indexing for faster search of RNA family databases. **BMC bioinformatics**, 9:131, Feb. 2008.

[103] T. Jiang, L. Wang, and K. Zhang. Alignment of trees - an alternative to tree edit. **Theor. Comput. Sci.**, 143(1):137–148, 1995.

[104] E. Y. Jin and C. M. Reidys. Asymptotic enumeration of RNA structures with pseudoknots. **Bulletin of mathematical biology**, 70:951–970, May 2008.

[105] T. Jörg, O. C. Martin, and A. Wagner. Neutral network sizes of biological RNA molecules can be computed and are not atypically small. **BMC bioinformatics**, 9:464, Oct. 2008.

[106] I. Kalvari, J. Argasinska, N. Quinones-Olvera, E. P. Nawrocki, E. Rivas, S. R. Eddy, A. Bateman, R. D. Finn, and A. I. Petrov. Rfam 13.0: shifting to a genome-centric resource for non-coding RNA families. **Nucleic acids research**, 46:D335–D342, Jan. 2018.

[107] R. Kleinkauf, M. Mann, and R. Backofen. antaRNA: ant colony-based RNA sequence design. **Bioinformatics (Oxford, England)**, 31:3114–3121, Oct. 2015.

[108] J. W. Klop, M. Bezem, and R. De Vrijer. Term rewriting systems. Cambridge University Press, 2001.

[109] R. V. Koodli, B. Keep, K. R. Coppess, F. Portela, E. participants, and R. Das. Eternabrain: Automated RNA design through move sets and strategies from an internet-scale RNA videogame. **PLOS Computational Biology**, 15(6):1–22, 06 2019.

[110] M. Kucharík, I. L. Hofacker, P. F. Stadler, and J. Qin. Basin hopping graph: a computational framework to characterize RNA folding landscapes. **Bioinformatics (Oxford, England)**, 30:2009–2017, July 2014.

[111] D. Lai, J. R. Proctor, and I. M. Meyer. On the importance of cotranscriptional RNA structure formation. **RNA**, 19(11):1461–1473, oct 2013.
[112] S. P. Lalley. Finite range random walk on free groups and homogeneous trees. **The Annals of Probability**, pages 2087–2130, 1993.

[113] K. Lari and S. J. Young. The estimation of stochastic context-free grammars using the inside-outside algorithm. **Computer Speech and Language**, 4:35–56, 1990.

[114] N. B. Leontis and E. Westhof. Geometric nomenclature and classification of RNA base pairs. **RNA**, 7 (4):499–512, 2001.

[115] A. Lescoute, N. B. Leontis, C. Massire, and E. Westhof. Recurrent structural RNA motifs, isostericity matrices and sequence alignments. **Nucleic acids research**, 33:2395–2409, 2005.

[116] 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–52, 2012.

[117] 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.

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

[119] R. Lorenz, C. Flamm, and I. L. Hofacker. 2d projections of RNA folding landscapes. In German Conference on Bioinformatics 2009, pages 11–20, Martin Luther University Halle-Wittenberg, Germany, Sept. 2009.

[120] W. A. Lorenz, P. Clote, and Y. Ponty. Asymptotics of RNA shapes. **Journal of Computational Biology**, 15(1):31–63, 2008.

[121] Z. J. Lu, J. W. Gloor, and D. H. Mathews. Improved RNA secondary structure prediction by maximizing expected pair accuracy. **RNA**, 15(10):1805–1813, Oct. 2009.

[122] 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.

[123] R. B. Lyngsø and C. N. Pedersen. RNA pseudoknot prediction in energy-based models. **Journal of computational biology : a journal of computational molecular cell biology**, 7:409–427, 2000.

[124] R. Lyngsø. Complexity of pseudoknot prediction in simple models. In Proceedings of ICALP, 2004.

[125] R. B. Lyngsø, J. W. Anderson, E. Sizikova, A. Badugu, T. Hyland, and J. Hein. FRNAkenstein: multiple target inverse RNA folding. **BMC Bioinformatics**, 13:260, 2012.

[126] J. Mañuch, C. Thachuk, L. Stacho, and A. Condon. NP-completeness of the energy barrier problem without pseudoknots and temporary arcs. **Natural Computing**, 10(1):391–405, nov 2010.

[127] D. H. Mathews. Using an RNA secondary structure partition function to determine confidence in base pairs predicted by free energy minimization. **RNA**, 10(8):1178–1190, jul 2004.

[128] D. H. Mathews and D. H. Turner. Dynalign: an algorithm for finding the secondary structure common to two RNA sequences. **Journal of Molecular Biology**, 317(2):191–203, mar 2002.

[129] D. H. Mathews, M. D. Disney, J. L. Childs, S. J. Schroeder, M. Zuker, and D. H. Turner. Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure. **Proc Natl Acad Sci U S A**, 101(19):7287–92, 2004.

[130] J. S. McCaskill. The equilibrium partition function and base pair binding probabilities for RNA secondary structure. **Biopolymers: Original Research on Biomolecules**, 29(6-7):1105–1119, 1990.

[131] Z. Miao, RNA-Puzzles Consortium, and E. Westhof. RNA-Puzzles Round III: 3D RNA structure prediction of five riboswitches and one ribozyme. **RNA**, 23:655–672, 2017.

[132] J. Michalik. *Non-redundant sampling in RNA bioinformatics*. PhD thesis, Interface Graduate school – Université Paris-Saclay, 2019.

[133] J. Michalik, 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.

[134] F. Michel and E. Westhof. Modelling of the three-dimensional architecture of group i catalytic introns based on comparative sequence analysis. **Journal of Molecular Biology**, 216(3):585–610, dec 1990.
[135] I. Miklós. *Computational Complexity of Counting and Sampling*. Discrete Mathematics and Its Applications. Chapman and Hall/CRC, feb 2019.

[136] I. Miklós, I. M. Meyer, and B. Nagy. Moments of the boltzmann distribution for RNA secondary structures. **Bull Math Biol**, 67(5):1031–1047, Sep 2005.

[137] M. Mishna and A. Rechnitzer. Two non-holonomic lattice walks in the quarter plane. **Theoretical Computer Science**, 410(38-40):3616–3630, 2009.

[138] M. Mohri. Semiring frameworks and algorithms for shortest-distance problems. **J. Autom. Lang. Comb.**, 7(3):321–350, Jan. 2002.

[139] M. E. Nebel. Combinatorial properties of RNA secondary structures. **Journal of Computational Biology**, 9(3):541–573, 2002.

[140] Y. Nesterov and A. Nemirovskii. *Interior-point polynomial algorithms in convex programming*, volume 13. Siam, 1994.

[141] M. E. Newman, S. H. Strogatz, and D. J. Watts. Random graphs with arbitrary degree distributions and their applications. **Phys. Rev. E**, 64(2):026118, August 2001.

[142] C. Nicaud and D. Gouyou-Beauchamps. Random generation using binomial approximations. **Discrete Mathematics & Theoretical Computer Science**, 2010.

[143] R. Nussinov and A. B. Jacobson. Fast algorithm for predicting the secondary structure of single-stranded RNA. **Proceedings of the National Academy of Sciences of the United States of America**, 77:6309–6313, Nov. 1980.

[144] C. Pivoteau. *Génération aléatoire de structures combinatoires : méthode de Boltzmann effective*. phdthesis, Université Pierre et Marie Curie - Paris VI, 2008.

[145] C. Pivoteau, B. Salvy, and M. Soria. Algorithms for combinatorial structures: Well-founded systems and newton iterations. **Journal of Combinatorial Theory, Series A**, 119(8):1711–1773, nov 2012.

[146] Y. Ponty. *Models for structured genomic sequences, random generation and applications*. Theses, Université Paris Sud - Paris XI, Nov. 2006.

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

[148] Y. Ponty. Efficient sampling of RNA secondary structures from the Boltzmann ensemble of low-energy: The boustrophedon method. **Journal of Mathematical Biology**, 56(1-2):107–127, 2008.

[149] Y. Ponty. Rule-weighted and terminal-weighted context-free grammars have identical expressivity. Research report, 2012.

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

[151] Y. Ponty, M. Termier, and A. Denise. GenRGenS: Software for Generating Random Genomic Sequences and Structures. **Bioinformatics**, 22(12):1534–1535, 2006.

[152] J. G. Propp and D. B. Wilson. Exact sampling with coupled markov chains and applications to statistical mechanics. **Random Structures & Algorithms**, 9(1-2):223–252, 1996.

[153] 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.

[154] J. Reeder and R. Giegerich. Design, implementation and evaluation of a practical pseudoknot folding algorithm based on thermodynamics. **BMC Bioinformatics**, 5:104, 2004.

[155] J. Reeder and R. Giegerich. Consensus shapes: an alternative to the sankoff algorithm for RNA consensus structure prediction. **Bioinformatics (Oxford, England)**, 21:3516–3523, Sept. 2005.

[156] C. M. Reidys, F. W. D. Huang, J. E. Andersen, R. C. Penner, P. F. Stadler, and M. E. Nebel. Topology and prediction of RNA pseudoknots. **Bioinformatics**, 27(8):1076–1085, Apr 2011.

[157] 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.
[158] 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.

[159] 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.

[160] 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.

[161] 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.

[162] K.-J. Räihä and E. Ukkonen. The shortest common supersequence problem over binary alphabet is NP-complete. *Theoretical Computer Science*, 16(2):187–198, 1981.

[163] 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.

[164] E. Rivas and S. Eddy. A dynamic programming algorithm for RNA structure prediction including pseudoknots. *J Mol Biol*, 285:2053–2068, 1999.

[165] C. Rovetta, J. Michálik, R. Lorenz, A. Tanzer, and Y. Ponty. Non-Redundant Sampling and Statistical Estimators for RNA Structural Properties at the Thermodynamic Equilibrium. working paper or preprint, Sept. 2019.

[166] F. Runge, D. Stoll, S. Falkner, and F. Hutter. Learning to design RNA. In International Conference on Learning Representations, 2019.

[167] D. Sankoff. Simultaneous solution of the RNA folding, alignment and protosequence problems. *SIAM Journal on Applied Mathematics*, 45(5):810–825, 1985.

[168] C. Saule, M. Régnier, J.-M. Steyaert, and A. Denise. Counting RNA pseudoknotted structures. *Journal of computational biology: a journal of computational molecular cell biology*, 18:1339–1351, Oct. 2011.

[169] G. Sauthoff, M. Möhl, S. Janssen, and R. Giegerich. Bellman's gap—a language and compiler for dynamic programming in sequence analysis. *Bioinformatics*, 29(5):551–560, 2013.

[170] A. Scala, L. Nunes Amaral, and M. Barthélémy. Small-world networks and the conformation space of a short lattice polymer chain. *Europhys. Lett.*, 55(4):594–600, 2001.

[171] S. Schirmer and R. Giegerich. Forest alignment with affine gaps and anchors, applied in RNA structure comparison. *Theor. Comput. Sci.*, 483:51–67, 2013.

[172] M. Schnall-Levin, L. Chindelevitch, and B. Berger. Inverting the Viterbi algorithm: an abstract framework for structure design. In *Machine Learning, Proceedings of the Twenty 9, 2008*, pages 904–911, 2008.

[173] 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.

[174] 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.

[175] 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.

[176] M. J. Smola, G. M. Rice, S. Busan, N. A. Siegfried, and K. M. Weeks. Selective 2'-hydroxyl acylation analyzed by primer extension and mutational profiling (SHAPE-MaP) for direct, versatile and accurate RNA structure analysis. *Nat Protoc*, 10:1643–1669, 2015.

[177] A. Spasic, S. M. Assmann, P. C. Bevilacqua, and D. H. Mathews. Modeling RNA secondary structure folding ensembles using SHAPE mapping data. *Nucleic Acids Res*, 46(1):314–323, 2017.
[178] J. Stombaugh, C. L. Zirbel, E. Westhof, and N. B. Leontis. Frequency and isostericity of RNA base pairs. **Nucleic acids research**, 37:2294–2312, Apr. 2009.

[179] D. Sundfeld, J. H. Havgaard, A. C. M. A. de Melo, and J. Gorodkin. Foldalign 2.5: multithreaded implementation for pairwise structural RNA alignment. **Bioinformatics**, 32(8):1238–1240, dec 2015.

[180] D. Surujon, Y. Ponty, and P. Clote. Small-world networks and RNA secondary structures. **Journal of computational biology: a journal of computational molecular cell biology**, 26(1):16–26, Jan. 2019.

[181] J. E. Tabaska, R. B. Cary, H. N. Gabow, and G. D. Stormo. An RNA folding method capable of identifying pseudoknots and base triples. **Bioinformatics**, 14(8):691–699, 1998.

[182] M. K. Takahashi and J. B. Lucks. A modular strategy for engineering orthogonal chimeric RNA transcription regulators. **Nucleic Acids Research**, 41(15):7577–7588, 2013.

[183] C. Theis, S. Janssen, and R. Giegerich. Prediction of RNA secondary structure including kissing hairpin motifs. In *Algorithms in Bioinformatics*, pages 52–64, Berlin, Heidelberg, 2010. Springer Berlin Heidelberg. ISBN 978-3-642-15294-8.

[184] I. Tinoco and C. Bustamante. How RNA folds. **J Mol Biol**, 293(2):271–81, Oct 1999.

[185] D. H. Turner and D. H. Mathews. NNDB: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure. **Nucleic acids research**, 38 part 2:(Database issue):D280–D282, Jan. 2010.

[186] V. Van Noort, B. Snel, and M. A. Huynen. The yeast coexpression network has a small-world, scale-free architecture and can be explained by a simple model. **EMBO Rep.**, 5(3):280–284, March 2004.

[187] G. Viennot and M. Vauchaussade de Chaumont. Enumeration of RNA secondary structures by complexity. In *Mathematics in Biology and Medicine*, pages 360–365, Berlin, Heidelberg, 1985. Springer Berlin Heidelberg. ISBN 978-3-642-93287-8.

[188] A. Viterbi. Error bounds for convolutional codes and an asymptotically optimum decoding algorithm. **IEEE Transactions on Information Theory**, 13(2):260–269, apr 1967.

[189] F. Voisin and M.-C. Gaudel. Drawing uniformly at random in dynamic sets of paths. HAL preprint, Oct. 2019.

[190] B. Voss, R. Giegerich, and M. Rehmsmeier. Complete probabilistic analysis of RNA shapes. **BMC biology**, 4:5, Feb. 2006.

[191] J. Waldispühl, S. Devadas, B. Berger, and P. Clote. Efficient algorithms for probing the RNA mutation landscape. **PLoS computational biology**, 4:e1000124, Aug. 2008.

[192] 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.

[193] 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.

[194] M. Waterman. Secondary structure of single-stranded nucleic acids. **Advances in Mathematics: Supplementary Studies**, 1:167–212, 1978.

[195] D. J. Watts and S. H. Strogatz. Collective dynamics of 'small-world' networks. **Nature**, 393(6684): 440–442, June 1998.

[196] S. E. Wells, P. E. Hillner, R. D. Vale, and A. B. Sachs. Circularization of mRNA by eukaryotic translation initiation factors. **Molecular cell**, 2:135–140, July 1998.

[197] H. S. Wilf. A unified setting for sequencing, ranking, and selection algorithms for combinatorial objects. **Advances in Mathematics**, 24:281–291, 1977.

[198] K. A. Wilkinson, E. J. Merino, and K. M. Weeks. Selective 2'-hydroxyl acylation analyzed by primer extension (SHAPE): quantitative RNA structure analysis at single nucleotide resolution. **Nat Protoc**, 1:1610–1616, 2006.

[199] S. Will, T. Joshi, I. L. Hofacker, P. F. Stadler, and R. Backofen. LocaRNA-p: accurate boundary prediction and improved detection of structural RNAs. **RNA (New York, N.Y.)**, 18:900–914, May 2012.
[200] S. Will, C. Otto, M. Miladi, M. Möhl, and R. Backofen. Sparse: quadratic time simultaneous alignment and folding of RNAs without sequence-based heuristics. **Bioinformatics (Oxford, England)**, 31: 2489–2496, Aug. 2015.

[201] M. T. Wolfinger, W. A. Svrcek-Seiler, C. Flamm, I. L. Hofacker, and P. F. Stadler. Efficient computation of RNA folding dynamics. **Journal of Physics A: Mathematical and General**, 37(17):4731–4741, apr 2004.

[202] A. R. Woods. Coloring rules for finite trees, and probabilities of monadic second order sentences. **Random Structures & Algorithms**, 10(4):453–485, 1997.

[203] S. Y. Wu, G. Lopez-Berestein, G. A. Calin, and A. K. Sood. RNAi therapies: Drugging the undruggable. **Science Translational Medicine**, 6(240):240ps7, 2014.

[204] S. Wuchty. Small worlds in RNA structures. **Nucleic Acids. Res.**, 31(3):1108–1117, February 2003.

[205] A. Xayaphoummine, T. Bucher, F. Thalmann, and H. Isambert. Prediction and statistics of pseudoknots in RNA structures using exactly clustered stochastic simulations. **Proc. Natl. Acad. Sci. U. S. A.**, 100 (26):15310–15315, 2003.

[206] 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.

[207] A. M. Yoffe, P. Prinsen, W. M. Gelbart, and A. Ben-Shaul. The ends of a large RNA molecule are necessarily close. **Nucleic acids research**, 39:292–299, Jan. 2011.

[208] J. N. Zadeh, B. R. Wolfe, and N. A. Pierce. Nucleic acid sequence design via efficient ensemble defect optimization. **Journal of Computational Chemistry**, 32(3):439–52, 2011.

[209] 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.

[210] 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.

[211] C. H. zu Siederdissen, S. J. Prohaska, and P. F. Stadler. Algebraic dynamic programming over general data structures. **BMC bioinformatics**, 16(19):S2, 2015.

[212] M. Zuker and D. Sankoff. RNA secondary structures and their prediction. **Bulletin of mathematical biology**, 46(4):591–621, 1984.

[213] M. Zuker and P. Stiegler. Optimal computer folding of large RNA sequences using thermodynamics and auxiliary information. **Nucleic Acids Research**, 9(1):133–148, 1981.