# Computational Complexity, Algorithmic Scope, and Evolution

Leonhard Sidl, Maximilian Faissner, Manuel Uhlir, Cristian A Velandia-Huerto, Maria Waldl, Hua-Ting Yao, Ivo L Hofacker, Peter F Stadler

► To cite this version:

Leonhard Sidl, Maximilian Faissner, Manuel Uhlir, Cristian A Velandia-Huerto, Maria Waldl, et al. 2025. hal-04715891v2

HAL Id: hal-04715891

https://hal.science/hal-04715891v2

Preprint submitted on 14 Feb 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Computational Complexity, Algorithmic Scope, and
Evolution

Leonhard Sidl<sup>1,2,‡</sup>, Maximilian Faissner<sup>1,‡</sup>, Manuel Uhlir<sup>1,‡</sup>, Cristian A. Velandia-Huerto<sup>1,3</sup>, Maria Waldl<sup>4</sup>, Hua-Ting Yao<sup>1</sup>, Ivo L. Hofacker<sup>1,2</sup>, Peter F. Stadler<sup>4,5,1,6§</sup>

E-mail: studla@bioinf.uni-leipzig.de

**Abstract.** Biological systems are widely regarded as performing computations. It is much less clear, however, what exactly is computed and how biological computation fits within the framework of standard computer science. Here we explore the idea that evolution confines biological computation to subsets of instances that can be solved efficiently with algorithms that are “hardcoded” in the system itself. We use RNA secondary structure prediction as a simple surrogate for developmental programs to demonstrate that the salient features of the genotype-phenotype map remain intact even if “simpler” algorithms are employed that correctly compute the structures only for small subsets of instances, albeit quantitative differences depending on the choice of alternative algorithms can be observed.

*Keywords:* Computational Complexity; Biological Computation; RNA Folding; Developmental Program; Genotype-Phenotype Map;

Submitted to: *J. Phys.: Complexity*

<footnote n="1">¹ Department of Theoretical Chemistry, University of Vienna, Währinger Straße 17, A-1090, Wien, Austria</footnote>

<footnote n="2">² Research Group Bioinformatics and Computational Biology, University of Vienna, Währinger Straße 29, A-1090, Wien, Austria</footnote>

<footnote n="3">³ Center for Anatomy and Cell Biology, Medical University of Vienna, Schwarzspanierstraße 17, 1090, Wien, Austria</footnote>

<footnote n="4">⁴ Bioinformatics Group, Department of Computer Science and Interdisciplinary Center for Bioinformatics, University Leipzig, D-04107 Leipzig, Germany</footnote>

<footnote n="5">⁵ Max Planck Institute for Mathematics in the Sciences, D-04103 Leipzig, Germany</footnote>

<footnote n="6">⁶ Santa Fe Institute, Santa Fe, NM 87501, USA</footnote>
# 1. Introduction

Biological systems are often perceived to perform computations [1] that solve complex problems at low computational cost as a consequence of evolutionary adaptation [2, 3]. In particular, the central “information metabolism” of a cell, i.e., DNA replication, transcription, and translation of RNA to proteins, is commonly described as a computation. Beyond these transformations of encoded information, and information processing in neural systems, it is often unclear what exactly a biological system is computing. Gene regulation has been described in terms of regulatory circuits akin to the hardware of a computing machine since the seminal work of Jacob and Monod [4, 5, 6]. Nevertheless, very little seems to be known about the computational problem(s) that this hardware is used to solve. Attempts to formalize regulatory mechanisms in terms of “I/O-maps” similar to logic gates [7] and the concept of a developmental program in evolutionary biology [8, 9, 10] have contributed to this perspective.

Despite the general agreement that “life computes”, the exact nature of natural computation is far from being well understood [11]. So far, biological computation has been studied predominantly from an information-theoretic and/or dynamical systems perspective [12, 13, 3]. Connections to theoretical computer science so far primarily concern the computational power of natural systems, such as chemical reaction networks [14], DNA computing [15, 16] or histone modification systems [17]. Typically, this line of work treats biological systems as paradigms for abstract models of computation, aiming to establish Turing universality or define computational limits. However, such general statements do not address how “difficult” the computational problems solved by biological systems actually are.

Here, we take a different point of view and explore whether there is a meaningful connection between the ability of a biological system to solve specific tasks and the concept of “complexity” in theoretical computer science, which quantifies the difficulty of computational problems and the resource requirements of the algorithms used to solve them. In a nutshell, a computational problem is an infinite set of individual instances, each with a defined solution, solved by algorithms that produce the appropriate answer from a given instance. The complexity of a problem is determined by the most difficult, i.e., most resource-consuming, instances. For example, the graph 3-coloring problem, a cornerstone of computational complexity studies, involves determining whether the vertices of a graph can be colored with three colors such that no two adjacent vertices share the same color. It is one of the classical NP-complete problems which are thought to require exponential running time [18, 19]. Restricting the set of instances, such as considering only triangle-free planar graphs, can drastically reduce its complexity, making it solvable in linear time [20].

Biological systems may likewise evolve to avoid inherently “hard” instances, focusing their solutions on specific subsets of problems. Levinthal’s paradox on protein folding observes that natural proteins fold remarkably fast into their native conformations, despite the immense size of the conformation space of peptide chains, which precludes exhaustive searching [21]. This led to the concept of folding funnels and energy landscapes that facilitate efficient folding [22, 23]. However, computational models often classify protein folding as an NP-complete problem [24, 25], see also [26]. Furthermore, many intrinsically disordered
proteins (IDPs) do not adopt a single stable three-dimensional structure [27].

It appears that nature does not solve the protein folding problem for all amino acid sequences. Instead, the problem is restricted, *by natural selection*, to a subset of amino acid sequences (instances) with suitable folding properties. We hypothesize that this may be a much more common phenomenon that pertains to biological information processing in general. The scope of computational problems in biology may be delineated by evolution in such a way that difficult instances are avoided. For example, given a developmental program, evolution may avoid genotypes (instances) that are “too difficult” for the machinery of development to correctly and efficiently translate into a viable phenotype. Such genotypes are unviable and are therefore selected against, effectively excluding these instances from the specification of the computational problem.

In this contribution we explore the impact of restricting the scope of a developmental program on evolution, using RNA secondary structure prediction as an example, since it has been widely used as a model system for simulating genotype-phenotype relationships [26, 29, 30, 31]. We introduce several different heuristics, each capturing different aspects of the folding process, and use the standard secondary structure model implemented by RNAfold in the ViennaRNA package [32] to define the ground truth of correctly computed structures. We shall see that the simplified algorithms solve only subsets of the instances, which in some cases are very small. Despite the drastic reduction of the accessible space of instances, we observe that the ability to evolve specific target structures remains intact. Additionally, only moderate quantitative effects on the evolutionary time-scales of reaching a target were found. The results presented here provide insight into how natural selection and computational constraints might jointly shape biological systems and lend credibility to the hypothesis that evolution can act to avoid computation of difficult instances.

## 2. Problem Complexity

In this section we briefly review the notion of problem complexity in theoretical computer science. Formally, a problem $\mathcal{P}$ is an infinite set of concrete realizations, called *instances*, together with a (possibly empty) set of solutions $S(\mathcal{J})$ for each instance $\mathcal{J} \in \mathcal{P}$. A particular instance $\mathcal{J}$ can thus be seen as the input for $\mathcal{P}$. For example, a specific RNA sequence $\mathcal{J}$ is an instance of the RNA folding problem, with the secondary structure $S(\mathcal{J})$ being the corresponding solution. The RNA FOLDING PROBLEM, thus, encompasses the computation of the secondary structure for any input sequence. The size of instance, $|\mathcal{J}|$, measures the amount of information necessary to specify it. In the example of RNA folding, $|\mathcal{J}|$ is simply the length of the input sequence $\mathcal{J}$. It is generally assumed that computing $S(\mathcal{J})$ from $\mathcal{J}$ will require more effort for larger instances.

An *algorithm* $\mathcal{A}$ that solves $\mathcal{P}$ takes any $\mathcal{J}$ as input and computes $S(\mathcal{J})$ using a certain amount of computational resources, which we denote by $\text{cost}(\mathcal{A}, \mathcal{J})$. We denote by $\mathfrak{A}[\mathcal{P}]$ the set of algorithms that solve $\mathcal{P}$, i.e., computes $S(\mathcal{J})$ for all $\mathcal{J} \in \mathcal{P}$. Algorithms, in turn, are specified with respect to some model of computation, such as the Turing machines conceptually underlying our digital computers. A plethora of different models of computation
have been studied, some of which, such as DNA computing [15, 16], are derived from biological processes. Many of these are computationally universal, in the sense that they can simulate Turing machines and thus compute everything that can be computed by a Turing machine [33, 34]. For the purpose of this contribution, the concrete choice of an underlying model of computation is not relevant, although the specification and comparison of algorithms requires that we use a fixed model.

The task most commonly encountered is to determine the complexity of a given problem $\mathbb{P}$. There are two ways of approaching this question: The most direct approach is to explicitly describe a specific algorithm $\mathcal{A}$, prove that it correctly solves $\mathbb{P}$, and to analyse $\mathcal{A}$ to determine the computational effort cost($\mathcal{A}, J$) required to run $\mathcal{A}$ for any instance $J$ of size $|J| = n$. Alternatively, one may proceed by showing that a problem $\mathbb{P}$ is at least as complex (or at most as complex) as another problem $\mathbb{P}^*$ for which complexity results are known. This route is taken in particular to establish that a problem is NP-complete [18].

In order to abstract from specific computational devices, programming languages, compilers, and the concrete encoding of particular instances $J$, one usually is interested in the scaling of an upper bound of the worst case computational cost as a function of the size $|J|$ of the input, i.e., a function $c_{\mathcal{A}}(n)$ that satisfies

$$c_{\mathcal{A}}(n) \geq k \max_{J:|J|=n} \text{cost}(\mathcal{A}, J) \qquad (1)$$

for all $n \geq n_0$, where $k > 0$ and $n_0 > 0$ are arbitrary constants. One says that $\mathcal{A}$ requires $\mathcal{O}(c(n))$ effort, or simply that $\mathcal{A}$ is in complexity class $\mathcal{O}(c(n))$. The complexity of a problem is then given by the minimal worst-case complexity of any algorithm that solves $\mathbb{P}$, i.e.,

$$c_{\mathcal{P}}(n) = \inf_{\mathcal{A} \in \mathfrak{A}[\mathbb{P}]} c_{\mathcal{A}}(n) \qquad (2)$$

Restriction of a problem to a subset of instances makes a problem easier in general, since $\mathbb{P}' \subseteq \mathbb{P}$ implies $c_{\mathcal{P}'}(n) \leq c_{\mathcal{P}}(n)$. As mentioned in the introduction, it is an active field of research in theoretical computer science to determine subclasses of difficult problems that admit efficient solutions, in particular in the area of combinatorial optimization problems on graphs.

Information processing in biological systems can be expected, at least in many cases, to be the result of a given algorithm $\mathcal{A}$ that is “hardcoded” e.g. in a gene-regulatory network [35]. In the case of RNA structures, the physical process of forming base pairs and extending and modifying helices may be thought of as the pertinent algorithm. The natural question then is, which instances are solvable by the given algorithm $\mathcal{A}$. This question has been addressed at least for some (classes of) algorithms in the past. Probably the best-understood case is the canonical greedy algorithm $\mathcal{G}$, which operates on constrained linear optimization problems. More precisely, the set $\mathbb{P}$ containing instances of the form $J = (X, w, \mathcal{F})$ where $X$ is a finite set, $\mathcal{F} \subseteq 2^X$ is a set of “allowed” subsets, and $w: X \to \mathbb{R}^+$ is a weight function. The greedy algorithm $\mathcal{G}$ sorts the set $X$ in decreasing order w.r.t. $w$ and initializes an empty solution set $B$. Then it tests, for all $x \in X$ in this order, whether $B \cup \{x\} \in \mathcal{F}$. If so, $x$ is added to $B$, otherwise $x$ is discarded. $\mathcal{G}$ thus attempts to find a subset $B \in \mathcal{F}$ that maximizes the score $\sum_{x \in B} w(x)$. In
fact, adding base pairs stepwisely such that the base pairs form a feasible secondary structure (see below) and the free energy decreases as much as possible in each step, constitutes a greedy algorithm for RNA folding.

Naturally, one asks under what conditions on $J = (X, w, \mathfrak{F})$ we can be sure that $\mathfrak{G}$ correctly solves the problem of maximizing $w$. A famous result states that this is the case for an arbitrary choice of weights if and only if the set system $(X, \mathfrak{F})$ is a certain generalization of a matroid [36]. On the other hand, if $X$ is the vertex set of a graph $G = (X, E)$ and $\mathfrak{F}$ are the independent sets of the graph $G$ (that is, the subsets of $G$ that contain no pair of adjacent vertices), then the correct solutions are that of the weighted independent vertex set problem. This optimization problem is well-known to be NP-hard [37] and in general is not solved correctly by $\mathfrak{G}$. Moreover, restricting the instance set e.g. to $P_5$-free graphs, it is possible to find polynomial-time algorithms [38], which are, however, entirely unrelated to $\mathfrak{G}$. It is not difficult to see that our example of RNA secondary structure prediction is not solved by the greedy algorithm in general.

### 3. RNA folding as a computational model of evolution

In the following paragraphs we describe the RNA folding problem in some more detail and use it as an example to address the question whether algorithms that correctly solve different sets of instances can have an impact on evolution. We use here the setting of developmental programs, i.e., algorithms $\mathcal{A}$ that take a genotype $J$ as an input and produce a phenotype as an output. For simplicity, we assume that only correctly constructed phenotypes $S(J)$ are viable. It might be more natural to ask whether an algorithm $\mathcal{A}$ computes a sufficiently close approximation to $S(J)$ instead of just distinguishing between correct and incorrect answers. Many hard optimization problems with corresponding NP-complete decision problems also do not admit efficient accurate approximations [39]. At least for an initial study into the interplay of problem complexity, algorithm capabilities, and their evolutionary consequences, we therefore restrict ourselves to the simplest case and ignore details such as the quality of an approximation. In more realistic models, we expect of course that fitness will depend on some measure of approximation quality instead of just being a yes/no answer.

Development is itself a complex problem, and there is at present no complete model of even the simplest developmental process that would link a genotype directly to a phenotype. Structure formation of biopolymers, however, has been used extensively as a convenient, computationally tractable toy model of genotype-phenotype (GP) maps. This is in particular the case for RNA secondary structure, for which the GP maps $J \mapsto S(J)$ for fixed $n = |J|$ have been studied extensively in the 1990s [28, 29]. We therefore (re)use the paradigm here as well.

For every RNA sequence $J$ of length $|J|$ over the alphabet $\{A, C, G, U\}$ of nucleotides we define the solution $S(J)$ as the minimum free energy secondary structure, i.e., the contact map, as predicted by RNAfold [32]. The evolutionary dynamics of an evolving RNA population is essentially determined by the underlying GP map.

More precisely, one observes qualitatively the same evolutionary dynamics dominated
by diffusion on neutral networks and rapid adaptation on the transitions between neutral networks, independent of particular choice of a fitness function that evaluates phenotypes only [40]: The GP map of RNA secondary structures exhibits a high degree of neutrality, admitting extensive, essentially connected, neutral networks of sequences that fold into a common structure. On the other hand, the neutral networks of any two abundant structures almost touch somewhere in sequence space, i.e., for any two structures $\phi_1$ and $\phi_2$ there are very similar sequences $I_1$ and $I_2$ with $S(I_1) = \phi_1$ and $S(I_2) = \phi_2$. As a consequence, populations show a diffusive motion on the neutral network, while mutations that fall outside the neutral network explore previously unseen structures at a constant rate [41, 42]. This gives rise to evolutionary trajectories in phenotype space that show punctuated equilibria, i.e., long periods of stagnation with rapid intermittent innovations [41, 31].

In this setting we can explore the effect of restricting folding to a subset of sequences that produce the correct structure $S(I)$ with simpler heuristic folding algorithms $\mathcal{A}(I)$ that work correctly only on a (small) subset $\mathbb{P}_{\mathcal{A}}$ of $\mathbb{P}$. To this end we introduce several simple heuristics that emphasize different aspects of RNA folding. We then ask how the evolution is influenced by the different scope of an algorithm $\mathcal{A}$ and by costs associated with the choice of $\mathcal{A}$.

In the next section we describe several variants of simple RNA folding algorithms $\mathcal{A}$. Each is designed to accept an arbitrary sequence as input. We define the minimum free energy structures computed by the folding program RNAfold [32], as the ground truth, i.e., as the correct solution $S(I)$ for each input sequence. The problem solved by an alternative folding algorithm is therefore

$$ \mathbb{P}_{\mathcal{A}} := \{I \mid \mathcal{A}(I) = S(I)\} \qquad (3) $$

That is, the subset of RNA sequences for which $\mathcal{A}$ produces the same structure as RNAfold. Sequences $I$ with $\mathcal{A}(I) \neq S(I)$ are treated as failures and considered as not viable. We then proceed to investigate to what extent this affects the structure of the genotype-phenotype map and eventually evolutionary processes.

## 4. Alternative RNA Folding Algorithms

All algorithms described in this section use, if applicable, the energy model implemented in ViennaRNA, which is also used by RNAfold, with the Turner2004 parameters [43] for evaluating secondary structures and enforce the constraint that hairpin loops must contain at least three unpaired bases. We investigate five simple folding algorithms that belong to two broad classes. All of them are inspired by the dynamical process of RNA secondary structure formation rather than exact solution by means of dynamic programming [44], which in contrast is quite far removed from the physics of the folding process. Tab. 1 summarizes the characteristics of each folding algorithm described below. Additional details on the effects of different folding algorithms on genotype-phenotype maps and structure diversity can be found in Supp. A.

The first group of algorithms is intended to approximate co-transcriptional folding. Base pairs are formed while the RNA sequence grows from the 5'-end towards the 3'-end, Fig. 1A. The folding process occurs in two stages. In the first stage, the model tries to pair newly
<table><caption>Table 1: Summary of folding algorithm characteristics.</caption><thead><tr><th>Folding algorithm</th><th>Co-transcription</th><th>Full length</th><th>Energy model</th></tr></thead><tbody><tr><td>RNAfold</td><td></td><td>✓</td><td>✓</td></tr><tr><td>Look behind fold</td><td>✓</td><td></td><td></td></tr><tr><td>Basic co-fold</td><td>✓</td><td></td><td></td></tr><tr><td>Best helix co-fold</td><td>✓</td><td></td><td>✓</td></tr><tr><td>Folding rule</td><td></td><td>✓</td><td>✓</td></tr><tr><td>Beam search</td><td></td><td>✓</td><td>✓</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Figure 1: Greedy-like approaches to RNA folding either simulate co-transcriptional folding (A) or operate on the complete sequence (B). In the first case, structures are formed while the sequence grows from the 5' to the 3' end and are subsequently kept in place. *Look Behind Fold* (A) closes a helix as soon as a seed comprising three base pairs become available. The RNA *Folding rule* (B) does a stepwise insertion of the most stable helix that does not conflict with previously inserted ones.</figcaption></figure>

transcribed base(s) to form seed helices, which are then used as reference in the next stage to complete the folding. We consider three variants of this approach:

* *Look behind fold* introduces one transcribed base at each step and pairs it with the closest complementary base if a size three seed helix can be formed. *Look behind fold* thus is a greedy-like approach that tends to form short helices as soon as possible. Once seed helices are identified, the algorithm attempts to extend the helices outwards by adding additional base pairs on each side. The sequence is also zipped inwards by identifying base pair stacks in between the seeds.

* *Basic co-transcriptional fold (Basic co-fold)* works in two separate stages. In **Stage 1**, the algorithm moves along the sequence from 5' to 3' and visits one base of the sequence after another. For each base, the algorithm searches for complementary positions among previously visited bases starting at the 5' end. When a basepair is found, the algorithm checks if it can be extended to a helix of length at least six. Once the first such helix is found, it is formed, the search procedure is stopped, and the algorithm advances to
the next position. Such a locally stable helix is then frozen in place and cannot be
replaced by an energetically more beneficial structure later on. After this first pass over
the sequence is complete, the intermediate structure only contains long helices that are
assumed to be the energetic backbone of the structure. In **Stage 2**, the procedure is
repeated with a minimal helix length of four. The algorithm moves from 5' to 3' again,
skipping already paired positions. For each unpaired base, the algorithm searches for
complimentary positions, starting from the 5', and checks if a helix of length at least four
can be extended from those two seed positions. In this second stage, special care is taken
to avoid the formation of pseudoknots.

*   **Best helix co-transcriptional folding (Best helix co-fold)** is a variation of the previous algorithm that can also be split into two separate stages. In **Stage 1**, the algorithm moves over the sequence from 5' to 3'. However, it does not process each base individually, but considers the next three new bases in each step. All base pairs which contain at least one of the three currently considered bases and only positions on the 5' side of those are then enumerated. Each of those seeds is then extended to the longest helix possible. The helix which has a length of at least six and the most beneficial effect on the energy of the complete structure is then chosen and formed. After this, the procedure is repeated for the next three bases. As in **Basic co-fold**, those helices are fixed and can not be replaced later in the algorithm. The energy is calculated using RNAfold with default parameters. In **Stage 2 Best helix co-fold** uses RNAfold to complete the structure and add smaller helices. Here, the helices inserted in the first pass are used as hard constraints and thus persist in final structure [45]. The main difference to the previous algorithm is that we reveal several new bases at the same time and then find the most stable helix using the new bases. This avoids traps formed by long helices that are globally unfavorable.

The second group of algorithms folds the entire input sequence and aims to efficiently identify a near-optimal secondary structure by heuristically exploring a subset of potential conformations, Fig. 1B. Starting from an open chain as the initial structure state, the folding progresses by moving to the next conformation that provides the greatest improvement in the free energy, based on simple rules such as adding a base pair or a helix. At each step, the algorithm selects the most energetically favorable addition without allowing the introduction of crossing base pairs (pseudoknots). Here, we consider two variants:

*   The *RNA Folding rule* was proposed in [46] as a greedy approximation to kinetic folding. At each step, the algorithm adds the next most favorable helix to the structure. Once a helix is added, the algorithm does not reconsider previous choices or allow for backtracking, and the process terminates when no further helix can be added without increasing the total folding energy.

*   *Beam search* is similar to *RNA Folding rule* but only one base pair is added at each step. In addition, it employs a beam search strategy to track suboptimal intermediate states. At each step, the *k*-best candidate structures are stored among all possible candidates generated from the previously stored *k*-best states. If *k* = 1, the method reduces to a greedy algorithm based on individual base pairs, similar to the *RNA Folding rule*. Since
GGGGUCUUCUACAGAUUAUACGCCUGCAACCUUGUGGGCCGUUCCUGAACACCU
. ((((((((.................((........((((((((...........)))))))))..)))))))))
(((((((((.................(........((((((((...........)))))))))..)))))))))
(((........................((((((((...........)))))))))..))))
(((((((((.................(........((((((((...........)))))))))..)))))))))
(((........................((((((((...........................)))))))))..))))
(((...)))..(((...........)))..(((...)))..(((...))))..

RNAfold

Beam search

Folding rule

Best helix co-fold

Basic co-fold

Look behind fold

<figure><img src="image_3.png" /><figcaption>(a) RNAfold</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Beam search</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) Folding rule</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d) Best helix co-fold</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(e) Basic co-fold</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(f) Look behind fold</figcaption></figure>

Figure 2: Predictions of the different algorithms for a randomly selected input sequence *I*. **Top.** Structures are displayed in dot-bracket (Vienna) notations: each base pair is shown as a matching pair of parentheses, unpaired nucleotides are shown as dots. In this case, none of the simplified algorithms correctly computes the solution *S*(*I*) defined by RNAfold. **Below.** Conventional secondary structure drawings.

a greedy approach is unlikely to yield globally optimal structures, retaining *k* suboptimal intermediates significantly increases the probability of identifying the most energetically favorable structure. However, increasing *k* improves solution quality at the cost of higher computational complexity, so selecting *k* involves a trade-off between accuracy and running time. In this work, the value of *k* is set to 50.

A closer inspection of the structures produced by the different algorithms quickly reveals a few general trends indicating that they focus on different parts of the shape space. Fig. 2 shows the predicted structures of a randomly chosen instance *I*. The baseline structure given
<figure><img src="image_9.png" /><figcaption>Figure 3: Cumulative distribution of Hamming distances between a random instance $I_0$ and the first sequence $I^*$ that folds correctly (i.e., matches that RNAfold prediction) encountered along an adaptive walks minimizing $d(\mathcal{A}(I), S(I))$. Data are for a sequence length of $n = 100$ and a sample size of 10 000 sequences for each algorithm.</figcaption></figure>

by RNAfold combines both long and short range interactions with a plethora of complex features such as bulges and interior loops. Contrasting that, the simpler alternative algorithms tend to lack such complex structures and are marked by longer and locally more stable helices. This is not surprising, since all except Beam Search utilize the paradigm of greedy algorithms and may easily get trapped in local minima when finding new helices. A particular example here is that Folding rule forms the helix at (1,47) instead of (2,55) because of the greedy choice which leads to a significantly different structure. Despite the similarities of Basic co-fold and Best helix co-fold, the predicted structures at times differ drastically, and in the example in Fig. 2.

*Density of Correctly Solvable Instances.* Since secondary structures are defined in terms of base pairs, two structures for the same sequence are naturally compared by counting the number of distinct base pairs, i.e., the cardinality of the symmetric difference of their sets of base pairs: $d(S, S') := |(S \cup S') \setminus (S \cap S')|$. We refer to $d(S, S')$ as the base pair distance.

For our purposes, it is particularly important to evaluate to what extent the alternative algorithms compute the correct structure. For the empirical results outlined in the next section, this is defined as producing the same structure as RNAfold. A sizable fraction of correct folds is observed only for three of the five algorithms: *Beam search*, *Folding rule*, and *Best helix co-fold*.

In order to obtain a more detailed picture of instances correctly solvable ($\mathcal{A}(I) = S(I)$) by each algorithm, we evaluated whether instances with correct solutions can be found close to any randomly chosen instance $I_0$. To answer this question we consider the
problem of minimizing $d(\mathcal{A}(I), S(I))$ using an adaptive walk starting at randomly chosen sequences $I_0$. The walk performs random mutations with non-increasing distance and stops at $d(\mathcal{A}(I), S(I)) = 0$ or no mutation can be made (see Supp. B.1 for more implementation details). We found that the base pair distance reaches 0 after a moderate number of optimization steps. The lengths of these adaptive walks, i.e., the number of accepted steps, serve as a convenient quantification for the ease of finding a solvable instance. The results are summarized in Fig. 3. A very small number of steps suffices when using *Beam search*, *Folding rule*, and *Best helix co-fold*. We therefore conclude that for these algorithms the set of solvable instances $\mathcal{P}_\mathcal{A}$ is large and more or less densely and uniformly distributed within $\mathbb{P}$.

The situation is different for *Basic co-fold* and *Look behind fold*. While it is unlikely to find instances with correct solutions for *Basic co-fold* at random, it is not very difficult to produce such instances by means of adaptive walk. The ease with which we can find a solvable instance close to a random $I_0$ suggests that $\mathcal{P}_\mathcal{A}$ is still large and approximately uniformly distributed, albeit much less dense than for the three algorithms discussed above. In contrast, finding solvable instances for *Look behind fold* is difficult, even impossible for a small fraction of instances (data not shown), and requires very large distances from $I_0$, substituting a substantial fraction of the initial sequence.

# 5. Neutral Networks

In order to investigate the impact of different algorithms and their ability to correctly solve the folding problem, we consider here properties of the genotype-phenotype map that have a key influence on evolvability. In Supp. A, we briefly review features that are independent of the correctness of the secondary structures. Since these turn out to be very similar for all algorithms considered here, we attribute differences in the simulations primarily to differences in likelihood of correct predictions.

Evolvability is closely linked to the extent of neutral networks, which can be quantified by the length distribution of neutral paths and the covering radio [29]. A *neutral path* starts at an instance $I_0$ that folds correctly. From a sequence $I$, a candidate $I'$ is generated by mutating an unpaired position in $S(I)$ or the two bases of base pair in $S(I)$ such that the Hamming distance from the starting point is non-decreasing, $d_H(I_0, I') \ge d_H(I_0, I)$. After a certain number of attempts, a strict increase of the distance is enforced. A candidate is accepted if $\mathcal{A}(I') = S(I') = S(I_0)$. The path terminates with a correctly folded sequence $I^*$ if no further step can be made. The sequences along the path are therefore correctly predicted by $\mathcal{A}$ to fold into the same structure as the starting sequence $I_0$. The Hamming distance $d_H(I_0, I^*)$ is an important characteristic of neutral networks. As shown in [29], such neutral paths can be expected to extend across the entire sequence space, i.e., $d_H(I_0, I^n) \approx n$ for RNAfold. For additional technical details on the analysis of neutral paths we refer to Supp. B.2.

Distributions of neutral path lengths for the different algorithms are summarized in Fig. 4a. We observe that they behave similar to RNAfold, indicating that in general there are still extensive neutral networks. Occasionally, we observe very short paths, although
<figure><img src="image_10.png" /><figcaption>Figure 4: Neutral paths as a function of the algorithm $\mathcal{A}$. (a) The distribution of maximum reached Hamming distances of a walk from a given start sequence. The distribution is drawn from 10 000 random start sequences of length 100. (b) The distribution of the residual Hamming distances of a walk for a given pair of start and target sequence. The distribution is obtained for 100 random target for each of 50 random start sequences of length 100.</figcaption></figure>

this appears with a non-negligible probability for *Beam search* and *Folding Rule*. For these, extending neutral paths is immediately infeasible for about 0.1% of the randomly sampled start sequences. In this case, none of the neighbors of $I_0$ correctly folds into the same structure. Interestingly, neutral paths almost never stops at low distances if initial steps are feasible.

A variation on the construction of neutral paths considers walks from $I_0$ towards a target sequence $I^*$ that is correctly computed by $\mathcal{A}$. In each step, the $d_H(I, I^*)$ is non-increasing and, after a certain number of steps, a strict decrease of the Hamming distance is enforced. Here, we are interested in the *residual Hamming distance* that remains when no further step towards the target is possible. Fig. 4b summarizes the final residual Hamming distances. More extensive neutral networks will approach the target more closely. This residual distance, however, is bounded by incompatibilities of $S(I_0)$ and the target sequence $I^*$, i.e., base pairs in $S(I_0)$ that cannot be formed by target sequence explain the sizable residual distance observed in Fig. 4b.

Even though for some algorithms only a small fraction of the sequences are correctly predicted, we observe very similar statistics for their neutral paths. The main difference, therefore, seems to be the difficulty of finding a viable, i.e., correctly predicted starting sequence $I_0$.

## 6. Simulated Evolutionary Trajectories

In order to evaluate the impact of the folding algorithms on an evolutionary process more directly we simulated a population of RNAs, initialized with random sequences that have correctly predicted structures, in which sequences replicate with structure-dependent fitness.
The fitness function is defined in terms of the similarity of $\mathcal{A}(I)$ with a pre-defined target structure $S^*$, provided the structure prediction is correct. Otherwise, we assume that $I$ is not viable, has fitness 0, and hence never enters the population. We use the fitness function

$$f(I) = \begin{cases} (0.01 + d(\mathcal{A}(I), S^*)/n)^{-1} & \text{if } \mathcal{A}(I) = S(I) \\ 0 & \text{if } \mathcal{A}(I) \neq S(I) \end{cases} \qquad (4)$$

in all simulations where $d$ denotes the base pair distance as above.

The simulations are performed in a stochastic flow reactor model based on the Gillespie algorithm [47, 48] and adapted to RNA evolution in [49, 31, 50]. A more detailed description can be found in Supp. B.3. We assume that only the replication event itself takes time, while mutations are introduced stochastically into the offspring as part of this process. Unviable offspring, i.e., sequences for which an algorithm $\mathcal{A}$ does not correctly predict the structure, are never introduced into the population. Such an unsuccessful replication event, however, takes time and causes the same exponentially distributed advance of the simulation clock as a successful event. In order to keep the population size constant, the offspring of successful replication event replaces a member of the population with uniform probability. The simulation terminates if half of the population folds into the target structure $S^*$.

<figure><img src="image_11.png" /><figcaption>Figure 5: Time required for population to reach average base pair distance of 30, 20, and 10 bp and the end of simulated evolutionary using two bacterial RNAs as targets: tRNA-Val(TAC) (left) and riboswitch aptamer DUF1646 (right). Data are taken from 50 simulation runs for each pair of target and folding algorithm. Values beyond 1.5 times of interquartile range (IQR) from the first or third quartile are considered as outliers (empty circles). *Look behind fold* failed to reach both targets, *Basic co-fold* could not reach the aptamer structure in any replicate within $10^7$ replication events, indicated by * and +, respectively.</figcaption></figure>

Fig. 5 shows the time in the evolution simulation at which a population reaches different average distances to the target structure, as well as the end time of the simulation, where 50% of the population has reached the target. Here, we consider two RNAs as targets: tRNA-Val(TAC) from *Escherichia coli* K-12 (76 nucleotides) and a Na+ riboswitch aptamer (DUF1646) from *Lactobacillus crispatus* ST1 (57 nucleotides). In general terms,
measuring the time required to reach a given base pair distance from the target yielded minor
differences among algorithms. However, *Basic co-fold* and *Look behind fold* exhibited slower
performance as the evolutionary trajectory approaches the target. For tRNA, *Look behind
fold* does not reach the target structure, for the ribozyme target, both *Basic co-fold* and *Look
behind fold* fail to find the target structure. Instead, the simulations get stuck at small residual
distances, suggesting that the target structure is either exceedingly rare or cannot be formed
at all.

**7. Discussion**

The scope of *A* in systems such as development influences the genotype-phenotype map, which might translate into possibly drastic effects on evolutionary trajectories [29, 31, 51]. We used several alternatives of RNA folding algorithms to investigate to what extent the scope of an algorithm, i.e., the set of instances that it solves correctly, has an impact on the dynamics of evolution. Most directly, the choice of *A* influences the density of the subset of sequence space that is accessible to evolutionary optimization.

Surprisingly, algorithms that correctly solve a noticeable fraction of the instances seem
to perform very similarly to RNAfold, the ground-truth. In computer simulations we observed
very similar neutral paths and very similar timescales for populations to approach a target
secondary structure. Even algorithms that do not solve most instances correctly, such as
*Beam search* and *Folding Rule*, still yield extensive neutral paths and thus large neutral
networks. Reduced density of the neutral networks does not seem to have much impact on the
evolutionary dynamics. We observe that the evolutionary process can confine exploration to
the subset $\mathbb{P}_{\mathcal{A}}$ without a significant reduction in the efficiency of the evolutionary optimization
process. Apparently, the ability to explore large parts of the fitness landscape along neutral
networks is sufficient to ensure evolvability. We suggest that the same effect pertains to
Levinthal's paradox: protein evolution is not significantly hampered by simply ignoring amino
acid sequences that do not fold efficiently. There does seem to be a noticeable effect if correct
solutions become sparse. In the case of *Basic co-fold* and *Look behind fold* we observed the
same qualitative behavior as with other algorithms, albeit at somewhat expanded time scales.

We have restricted our attention here to the simplest scenario, considering only the scope
of the algorithms. We also restricted ourselves to “external” measures, comparing independent
evolution experiments. In a more elaborate setting, one could consider the direct competition
of two algorithms by associating with each individual in the population also a fixed algorithm
that is inherited in each replication event. In this setting it becomes possible to evaluate
 (relative) fitness contributions associated with the differences in scope of two algorithms.

At first glance, the fact that some algorithms very rarely solve a random instance correctly
might make them non-starters. However, we have used here a very strict setting in which
computational results $\mathcal{A}(I)$ need to be perfect solutions $S(I)$, declaring all other instances as
unviable. This may not be the case and approximations to the correct solution may already
yield some benefit, i.e., a non-zero fitness. In this case we can expect natural selection to
zoom in on nearby “easy” instances, which seem to be true according to Fig. 3.
We also have not included any direct fitness contributions associated with computation, such as research consumption for the prediction. While this may be less important for developmental programs, we suspect that the cost of computation becomes an important issue in neural systems, e.g., in the processing of external signals. What exactly are the trade-offs between generality of the problem sets and the cost of execution in this case? After all, we expect that in general an algorithm $\mathcal{A}'$ that solves a proper superset of $P[\mathcal{A}]$ will also require more resources. Our simulation data indicate that large fractions of instances that are not solvable by $\mathcal{A}$ do not hamper evolvability. This suggests that in general there is not much of a selective pressure to enlarge the scope of $\mathcal{A}$ imposed by the distribution of correct instances. In particular, significant fraction of unviable offsprings has little impact on evolvability and on the speed of evolutionary adaptation as long the generation of unviable offsprings does not incur a direct fitness cost for the parent.

In the case of protein folding, folding funnels may serve as an indication that there is selection against hard instances. Similarly, we observed in simple computer experiment (detailed in Supp. C) that RNA sequences that evolved under strong selection for a conserved structure are much more often predicted correctly by *Folding Rule*, the algorithm that was explicitly designed to capture the folding process. Similar effects were not observed for other algorithms that are arguably less close to the physics of RNA folding. While certainly not a proof, we interpret this result as an indication that the algorithmic scope may have a tangible impact on evolution.

From a biological point of view, moreover, we will eventually have to return to the question how one can determine the pertinent algorithms in detail, and how their resource consumption and input dependence could be quantified. It would appear that simplified computational models, including but not limited to the RNA folding models used here, can provide valuable insights, even though they probably fall short of being particularly realistic models of computation in biological systems.

Returning to the question whether the concepts developed in theoretical computer science to describe computational complexity can be useful in describing and understanding computation in biological systems, we arrive at an affirmative answer. However, instead of considering problems $P$ that are defined *a priori*, the focus will have to be on the algorithms $\mathcal{A}$ that have proved sufficient to ensure survival in evolutionary history. In order to make progress towards understanding *what* a biological system computes [52], it will be necessary to start from the *how*. Instead of problem complexity, the scope of “hardwired” algorithms is taking center stage. The formal framework of complexity theory, thus, remains an interesting approach, even though the types of questions are likely going to be different.

## Data Availability

The implementation of folding approaches and simulations can be found at https://github.com/ViennaRNA/VIECPLX.
Acknowledgements

Research in the Stadler lab is supported by the German Federal Ministry of Education and Research BMBF through DAAD project 57616814 (SECAI, School of Embedded Composite AI). PFS gratefully acknowledges discussions with Gülce Kardeş and David Wolpert in spring 2023 and feedback received on some of these ideas during the CSH Workshop *Trade-offs between thermodynamic cost, intelligence and fitness in living organisms*. HTY is funded by Austrian Science Fund (FWF), grant no. I 4520. MF, MU and MW are funded by the Novo Nordisk Foundation (grant NNF21OC0066551 “MATOMIC”). CAVH was funded by FWF grant number SFB F80-01. LS was funded by the FWF project number I 6440-N.

References

[1] Melanie Mitchell. Ubiquity symposium: Biological computation. *Ubiquity*, 2011:3, 2011.

[2] Christopher P Kempes, David Wolpert, Zachary Cohen, and Juan Pérez-Mercader. The thermodynamic efficiency of computations made in cells across the range of life. *Philos Trans A Math Phys Eng Sci*, 375(2109):20160343, 2017.

[3] David H. Wolpert, Chris Kempes, Peter F. Stadler, and Joshua A. Grochow, editors. *The Energetics of Computing in Life and Machines*. SFI Press, Santa Fe, NM, 2019.

[4] François Jacob and Jacques Monod. Genetic regulatory mechanisms in the synthesis of proteins. *J. Mol. Biol.*, 3(3):318–356, 1961.

[5] Douglas H Erwin and Eric H Davidson. The evolution of hierarchical gene regulatory networks. *Nat Rev Genet*, 10(2):141–148, 2009.

[6] Harold D. Kim, Tal Shay, Erin K. O’Shea, and Aviv Regev. Transcriptional regulatory circuits: Predicting numbers from alphabets. *Science*, 325:429–432, 2009.

[7] David C. Krakauer, Lydia Müller, Sonja J. Prohaska, and Peter F. Stadler. Design specifications for cellular regulation. *Th. Biosci.*, 135:231–240, 2016.

[8] G Oster and P Alberch. Evolution and bifurcation of developmental programs. *Evolution*, 36(3):444–459, 1982.

[9] Detlev Arendt, Jacob M. Musser, Clare V. H. Baker, Aviv Bergman, Connie Cepko, Douglas H. Erwin, Mihaela Pavlicev, Gerhard Schlosser, Stefanie Widder, Manfred D. Laubichler, and Günter P. Wagner. The origin and evolution of cell types. *Nat Rev Genet*, 17:744–757, 2016.

[10] Somya Mani and Tsvi Tlusty. A topological look into the evolution of developmental programs. *Biophys. J.*, 120(19):4193–4201, 2021.

[11] Dominique Chu, Mikhail Prokopenko, and J. Christian J. Ray. Computation by natural systems. *Interface Focus*, 8(6):20180058, 2018.

[12] Keyan Zahedi and Nihat Ay. Quantifying morphological computation. *Entropy*, 15:1887–1915, 2013.

[13] Neil Dalchau, Gregory Szép, Rosa Hernansaiz-Ballesteros, Chris P. Barnes, Luca Cardelli, Andrew Phillips, and Attila Csikász-Nagy. Computing with biological switches and clocks. *Natural Computing*, 17(4):761–779, 2018.

[14] François Fages, Guillaume Le Guludec, Olivier Bournez, and Amaury Pouly. Strong Turing completeness of continuous chemical reaction networks and compilation of mixed analog-digital programs. In J. Feret and H. Koepp, editors, *Computational Methods in Systems Biology. CMSB 2017*, volume 10545 of *Lect. Notes Comp. Sci.*, pages 108–127, Cham, 2017. Springer.

[15] David Soloveichik, Georg Seelig, and Erik Winfree. DNA as a universal substrate for chemical kinetics. *Proc. Natl. Acad. Sci. USA*, 107:5393–5398, 2010.

[16] Lulu Qian, David Soloveichik, and Erik Winfree. Efficient Turing-universal computation with DNA polymer. In Yasubumi Sakakibara and Yongli Mi, editors, *DNA Computing and Molecular*
*Programming*, volume 6518 of *Lect. Notes Comp. Sci.*, pages 123–140, Berlin, Heidelberg, 2011. Springer.

[17] B. Bryant. Chromatin computation. *PLoS ONE*, 7:e35703, 2012.

[18] Richard M. Karp. Reducibility among combinatorial problems. In R. E. Miller, J. W. Thatcher, and J. D. Bohlinger, editors, *Complexity of Computer Computations*, pages 85–103. Plenum, New York, 1972.

[19] Russell Impagliazzo and Ramamohan Paturi. The complexity of k-SAT. In *Proc. 14th IEEE Conf. on Computational Complexity*, pages 237–240, 1999.

[20] Zdeněk Dvořák, Ken-Ichi Kawarabayashi, and Robin Thomas. Three-coloring triangle-free planar graphs in linear time. *ACM Trans. Algorithms*, 7(4):41, 2011.

[21] Cyrus Levinthal. Are there pathways for protein folding? *J. Chim. Phys.*, 65:44–45, 1968.

[22] R Zwanzig, A Szabo, and B. Bagchi. Levinthal’s paradox. *Proc. Natl. Acad. Sci. USA*, 89:20–22, 1992.

[23] Martin Karplus. The Levinthal paradox: yesterday and today. *Folding and Design*, 2:S69–S75, 1997.

[24] R Unger and J Moult. Finding the lowest free energy conformation of a protein is an np-hard problem: proof and implications. *Bull Math Biol*, 55:1183–1198, 1993.

[25] B Berger and T. Leighton. Protein folding in the hydrophobic-hydrophilic (HP) model is NP-complete. *J Comput Biol*, 5:27–40, 1998.

[26] Christophe Guyeux, Nathalie M-L Côté, Jacques M Bahi, and Wojciech Bienia. Is protein folding problem really a NP-complete one? First investigations. *J Bioinform Comput Biol*, 12:1350017, 2014.

[27] Rakesh Trivedi and Hampapathalu Adimurthy Nagarajaram. Intrinsically disordered proteins: An overview. *Int J Mol Sci.*, 23(22):14050, 2022.

[28] Walter Fontana, Peter F Stadler, Erich G Bornberg-Bauer, Thomas Griesmacher, Ivo L. Hofacker, and Thomas Griesmacher, Manfred Ivo L Tacker, Pedro Tarazona, Edward D Weinberger, and Peter Schuster. RNA folding landscapes and combinatory landscapes. *Phys. Rev. E*, 47:2083–2099, 1993.

[29] Peter Schuster, Walter Fontana, Peter F Stadler, and Ivo L Hofacker. From sequences to shapes and back: A case study in RNA secondary structures. *Proc. Roy. Soc. Lond. B*, 255:279–284, 1994.

[30] Walter Fontana and Peter Schuster. Shaping space. The possible and the attainable in RNA genotype-phenotype mapping. *J. Theor. Biol.*, 194:491–515, 1998.

[31] Walter Fontana and Peter Schuster. Continuity in evolution. On the nature of transitions. *Science*, 280:1451–1455, 1998.

[32] Ronny Lorenz, Stephan H Bernhart, Christian Höner zu Siederdissen, Hakim Tafer, Christoph Flamm, Peter F. Stadler, and Ivo L. Hofacker. ViennaRNA Package 2.0. *Alg. Mol. Biol.*, 6:26, 2011.

[33] Udi Boker and Nachum Dershowitz. How to compare the power of computational models. In S. Barry Cooper, Benedikt Löwe, and Leen Torenvliet, editors, *New Computational Paradigms*, volume 3526 of *Lect. Notes Comp. Sci.*, pages 54–64, Berlin, Heidelberg, 2005. Springer.

[34] Naresh R. Shanbhag, Subhasish Mitra, Gustavode de Veciana, Michael Orshansky, Radu Marculescu, Jaijeet Roychowdhury, Douglas Jones, and Jan M. Rabaey. The search for alternative computational paradigms. *IEEE Design and Test*, 25:334–343, 2008.

[35] Eric H. Davidson. *The Regulatory Genome: Gene Regulatory Networks In Development And Evolution*. Academic Press, Burlington, MA, 2006.

[36] Paul Helman, Bernard M. E. Moret, and Henry D. Shapiro. An exact characterization of greedy structures. *SIAM J. Discrete Math.*, 6, 1993.

[37] M. R. Garey and D. S. Johnson. “Strong” NP-completeness results: Motivation, examples, and implications. *J. ACM*, 25(3):499–508, 1978.

[38] Daniel Lokshantov, Martin Vatshelle, and Yngve Villanger. Independent set in $P_5$-free graphs in polynomial time. In Chandra Chekuri, editor, *SODA ’14: Proceedings of the twenty-fifth annual ACM-SIAM symposium on Discrete algorithms*, pages 570–581, Philadelphia, PA, 2014. SIAM.

[39] David Zuckerman. On unapproximable versions of NP-complete problems. *SIAM J. Computing*, 25:1293–1304, 1996.

[40] Peter F. Stadler. Fitness landscapes arising from the sequence-structure maps of biopolymers. *J. Mol. Struct. (THEOCHEM)*, 463:7–19, 1999.

[41] Martijn A. Huynen, Peter F. Stadler, and Walter Fontana. Smoothness within ruggedness. The role of
neutrality in adaptation. *Proc. Natl. Acad. Sci. USA*, 93(1):397-401, 1996.

[42] Martijn A. Huynen. Exploring phenotype space through neutral evolution. *J. Mol. Evol.*, 43(3):165-169, 1996.

[43] David H Mathews, Matthew D Disney, Jessica L Childs, Susan J Schroeder, Michael Zuker, and Douglas H Turner. Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure. *Proc. Natl. Acad. Sci. USA*, 101(19):7287-7292, 2004.

[44] Michael Zuker and Patrick Stiegler. Optimal computer folding of larger RNA sequences using thermodynamics and auxiliary information. *Nucleic Acids Research*, 9(1):133-148, 1981.

[45] Ronny Lorenz, Ivo L. Hofacker, and Peter F. Stadler. RNA folding with hard and soft constraints. *Alg. Mol. Biol.*, 11:8, 2016.

[46] Hugo M Martinez. An rna folding rule. *Nucleic Acids Res*, 12:323-334, 1984.

[47] Daniel T Gillespie. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. *J. Comp. Phys.*, 22(4):403-434, 1976.

[48] Daniel T Gillespie. Exact stochastic simulation of coupled chemical reactions. *J. Phys. Chem.*, 81(25):2340-2361, 1977.

[49] Walter Fontana and Peter Schuster. A computer model of evolutionary optimization. *Biophysical chemistry*, 26(2-3):123-147, 1987.

[50] Andreas Wernitznig. RNA Optimization in Flow Reactors: A study in silico. PhD thesis, Univ. Vienna, 2001.

[51] Jacobo Aguirre, Javier M. Buldú, Michael Stich, and Susanna C. Manrubia. Topological structure of the space of phenotypes: The case of RNA neutral networks. *PLoS ONE*, 6(10):e26324, 2011.

[52] Sonja J. Prohaska, Peter F. Stadler, and Manfred Laubichler. How and what does a biological system compute? In David H Wolpert, Chris Kempes, Peter F. Stadler, and Joshua A. Grochow, editors, *The Energetics of Computing in Life & Machines*, pages 169-189. SFI Press, Santa Fe, 2019.