# mRNA secondary structure prediction using utility-scale quantum computers

Dimitris Alevras1, Mihir Metkar3, Takahiro Yamamoto2, Vaibhaw Kumar1, Triet
Friedhoff1,
Jae-Eun Park1, Mitsuharu Takeori2, Mariana LaDue1, Wade Davis3, and Alexey
Galda3  1IBM Quantum, New York, USA  2IBM Quantum, Tokyo, Japan  3Moderna,
Cambridge, USA

###### Abstract

Recent advancements in quantum computing have opened new avenues for tackling
long-standing complex combinatorial optimization problems that are intractable
for classical computers. Predicting secondary structure of mRNA is one such
notoriously difficult problem that can benefit from the ever-increasing
maturity of quantum computing technology. Accurate prediction of mRNA
secondary structure is critical in designing RNA-based therapeutics as it
dictates various steps of an mRNA life cycle, including transcription,
translation, and decay. The current generation of quantum computers have
reached utility-scale, allowing us to explore relatively large problem sizes.
In this paper, we examine the feasibility of solving mRNA secondary structures
on a quantum computer with sequence length up to 60 nucleotides representing
problems in the qubit range of 10 to 80. We use Conditional Value at Risk
(CVaR)-based VQE algorithm to solve the optimization problems, originating
from the mRNA structure prediction problem, on the IBM Eagle and Heron quantum
processors. To our encouragement, even with “minimal” error mitigation and
fixed-depth circuits, our hardware runs yield accurate predictions of minimum
free energy (MFE) structures that match the results of the classical solver
CPLEX. Our results provide sufficient evidence for the viability of solving
mRNA structure prediction problems on a quantum computer and motivate
continued research in this direction.

##  I Introduction

Messenger RNA (mRNA) plays a pivotal role in the central dogma of molecular
biology, acting as the intermediary between the genetic blueprint in DNA and
the functional molecule, protein. The functionality of mRNA is inherently
linked to its secondary structure, which arises from intramolecular base
pairing and folding patterns [1]. Understanding these structures is crucial
for elucidating the mechanisms of gene regulation, translation, and
degradation kinetics, and inter-biomolecular interactions [2]. Additionally,
the emergence of mRNA-based therapeutics has ignited a paradigm shift in the
landscape of medicine, offering novel solutions for the treatment of various
diseases [3]. This transition from understanding mRNA’s fundamental biological
role to harnessing its therapeutic potential marks a significant milestone in
biomedical research and holds immense promise for revolutionizing healthcare
practices. However, predicting mRNA secondary structures is a formidable
combinatorial problem due to the vast number of possible configurations and
their associated energy states [4]. Solving this problem accurately and
efficiently is essential for advancing our knowledge in genetics, aiding in
the design of mRNA-based drugs, and developing gene therapies.

The problem of RNA secondary structure prediction with pseudoknots is known to
be NP-complete [5]. State-of-the-art classical approaches for the secondary
structure prediction of mRNA have primarily relied on dynamic programming
algorithms, such as the Zuker algorithm [6] and its widely used
implementations, MFold [7] and ViennaRNA [8], to name a few. These methods
calculate the minimum free energy (MFE) structures by systematically
evaluating all possible base pairings and applying thermodynamic models to
estimate their stability [9]. Subsequent enhancements have incorporated
stochastic context-free grammars (SCFGs), as seen in tools like Infernal [10],
which align sequences to RNA structure profiles to predict conserved
structures across related sequences. The integration of machine learning
techniques was implemented in algorithms like ContextFold [11] that utilize
both local and global contextual features to improve prediction accuracy.
Despite these advancements, classical approaches are still constrained by the
combinatorial explosion of possible structures, lack of generalization to
novel RNAs [12], and the approximations required to make the problem tractable
[13]. As a result, there is a growing interest in leveraging quantum computing
to overcome these limitations and achieve more accurate predictions of mRNA
secondary structures.

Classical hardness of mRNA secondary structure prediction naturally positions
it as a compelling application for quantum computers. Despite still being
relatively scarce, recent developments have shown promising advancements in
this area. In Refs. [14, 15], quantum annealing-based hardware was used to
minimize the energy function that corresponds to the folding energy landscape
of mRNA molecules. However, quantum annealers are not universal quantum
processors, which limits the set of algorithms amenable to this technology
compared to the universal gate-based quantum computers. In another recent
paper [16], the Quantum Approximate Optimization Algorithm (QAOA) [17] was
used to solve the secondary structure prediction problem for toy-model RNA
sequences using 12 qubits on a simulator and 4 qubits on a gate-based quantum
processor.

Algorithms such as QAOA and Variational Quantum Eigensolver (VQE) fall under
the class of variational quantum algorithms that can be implemented on the
current generation of gate-based quantum computers. Guided by a classical
optimization feedback loop, the goal of such algorithms is to find the quantum
state that exhibits sufficient overlap with the ground state of the
optimization problem. Success of such schemes depends heavily on the choice of
a suitable ansatz, objective function, and efficient classical optimizer, as
highlighted in the literature. However, most of such studies are either based
on numerical simulations or examine relatively small problems on quantum
hardware, with only a few notable exceptions [18, 19, 20, 21, 22].

As universal quantum processors surpassed one hundred qubits, we have entered
the “utility” era of quantum computing [23], where quantum computers can
perform reliable computations on problem sizes beyond the range that can be
handled by brute-force classical schemes. Encouraged by these developments, in
this work we examine the feasibility of solving relatively longer mRNA
sequences on a universal utility-scale quantum processor. The goal here is to
demonstrate that variational quantum algorithms can be fully and successfully
executed on such hardware, with hundreds of quantum circuit calls during the
hyperparameter training step along with basic error mitigation, yielding
accurate results that match those obtained with classical optimization
schemes.

In this work, we use the Conditional Value at Risk Variational Quantum
Eigensolver (CVaR-VQE), proposed as a modification to the traditional VQE
algorithm, that employs CVaR values of the energy eigenfunctions as the
objective function, to achieve better convergence [24]. Additionally, we
choose Nakanishi-Fujii-Todo (NFT) [25] as a classical optimizer along with the
hardware-efficient “two-local” ansatz. The combination of the NFT optimizer
and two-local ansatz has been shown to be relatively robust against hardware
noise and exhibits better convergence properties [26].

First, we formulate the secondary structure prediction problem as a binary
optimization problem. We then use CPLEX [27] to solve the formulated problem.
We provide numerical evidence that as the problem size increases, classical
schemes like CPLEX show an extremely unfavorable runtime scaling. Next, we run
these problems using both quantum hardware and simulations, and present
results on RNA sequence lengths ranging from 15 to 42 bases representing
optimization problems in the range of 10 to 80 qubits. To the best of our
knowledge, our study is the first to examine RNA sequences up to 42
nucleotides on a quantum processor.

The paper is organized as follows. In Section II, we offer a formal definition
of the problem and give a brief description of the structures that provide
more or less stability. In Section III, we briefly describe the current
classical methods and describe state-of-the-art classical algorithms, with
focus on dynamic and mathematical programming approaches. We also provide our
formulation of the problem as a mathematical programming problem. Section IV
focuses on the quantum formulation of the problem and describes the CVaR-VQE
method used in this work. Section V contains the results of our experiments,
using both classical simulators and quantum processors. Finally, in Section
VI, we conclude with a summary of our results and an outlook.

Figure 1: An RNA secondary structure illustrating five kinds of secondary
structural elements: hairpin loops (blue), bulge (yellow), helix (black),
internal loop (green) and multi-branched loop (cyan) and a pseudoknot (pink)
[28].

##  II Problem definition

Figure 2: Example of an RNA sequence and its predicted MFE structure
demonstrating its folded state. Optimal folding calculated by ViennaRNA
RNAfold server. Graphic generated with ViennaRNA forna server [29].

RNA commonly is a polymer of four nucleotides (bases): adenine (A), cytosine
(C), guanine(G), and uracil (U). Under physiological conditions, RNA adopts
compact structures driven by hydrophobic, base-stacking, and base-pairing
interactions [1]. Typically these interactions facilitate the formation of an
A-form helix containing canonical Watson-Crick-Franklin base-pairs. The folded
RNA conformation can be depicted as secondary structures, comprising various
secondary structural elements that arise from base-pairing and base-stacking
interactions (see Fig. 1). Understanding RNA secondary structure is crucial
for gaining insights into its function and designing RNA-based therapeutics
[2]. The RNA secondary structure prediction problem is the problem of finding
the most stable folding of the sequence of bases [30]. The exact
characteristics of this “optimal” folding are not completely understood,
making approximate models well-suited for this purpose. In fact, in certain
cases, there may exist multiple optimal structures, underscoring the value of
methods that can produce numerous high-quality solutions.

According to thermodynamic principles, formation of stem-loop structure is
feasible if the entropic costs of closing loops and ordering nucleotides in a
helix are offset by free energy released due to base-stacking and pairing
interactions [31]. As a result, longer helices, or stems, are typically more
stabilizing for secondary structures, whereas loops and single-stranded
regions tend to introduce instability. To predict the most favorable folding,
a computational model must determine the structure with the lowest free energy
(MFE) compared to the unfolded state, which is the conformation that RNA is
most likely to adopt under physiological conditions [31].

##  III Classical approaches

The most common way of predicting secondary structure is to minimize the free
energy change upon moving from an unfolded state to the most thermodynamically
stable folded state, see Fig. 2 [32]. This secondary structure which has the
lowest possible energy is defined as the MFE structure. Because MFE reflects
the most thermodynamically stable structure, many approaches for predicting
secondary structure concentrate on identifying structures through MFE
minimization. However, this problem, viewed as an optimization problem, is
generally NP-hard. There are approximation schemes designed to optimize
simplified models by excluding specific structures like pseudoknots (see Fig.
1) and focusing only on a subset of possible structures [30, 33]. These
include dynamic programming algorithms (described below), along with genetic
algorithms [34] and machine learning algorithms [35], to name a few. See Ref.
[36] for a review of classical methods for RNA structure prediction.

The potential of quantum computing to tackle this problem stems from two key
aspects:

  * •

A Quadratic Unconstrained Binary Optimization (QUBO) problem can be converted
to the problem of finding the lowest energy eigenstate of an Ising Hamiltonian
and thus can be solved by a quantum algorithm;

  * •

The secondary structure prediction problem can be formulated as a QUBO.

Previous efforts have focused on developing an objective function to minimize
over binary variables, with subsequent optimization carried out on a quantum
annealing hardware [14, 15]. In this work, we use mathematical programming to
represent the problem in binary variables and then convert the resulting
program to QUBO, which we solve on a universal quantum computer.

###  III-A Dynamic programming algorithms

One approach to approximate the “optimal” folding is to maximize the number of
paired bases. This can be accomplished by the dynamic programming algorithm of
Nussinov [30]. The algorithm is based on the fact that the maximum number of
paired bases for sequence $i,\ldots,j$ can be calculated from the maximum
number of pairs for smaller length sequences. There are only four possible
ways to construct nested base-pairs from smaller sub-sequences [37]. It is
therefore possible to recursively calculate the maximum number of pairs. Even
though maximizing the number of paired bases is not sufficient, Nussinov’s
algorithm provides the basic mechanism of subsequent dynamic programming
algorithms, such as Zuker’s algorithm [38]; one has to replace the number of
maximum pairs in Nussinov’s algorithm with a scoring function (such as MFE)
and calculate the score for each sub-sequence. The best approximation to what
“optimal” means is the one that results in the best score (e.g., MFE) for the
sequence.

Dynamic programming algorithms do not account for all structures, especially a
specific type of structure known as a pseudoknot, shown in Fig. 1. While not
present in all mRNA sequences, and particularly relatively short sequences
considered in this work, pseudoknots are crucial for the function of several
important RNA elements, including regulation of translation and splicing [39]
and the binding of small molecules [40, 41, 42]. However, for many practical
purposes, approximation methods that do not explicitly account for pseudoknots
can still provide reasonable and efficient predictions of mRNA secondary
structures. The current classical algorithms utilized in software achieve
65-73% accuracy in predicting secondary structures observed in the laboratory
[43].

###  III-B Mathematical programming models

Among classical schemes, mathematical programming models and integer linear
programming, in particular, are often used to solve the RNA folding problem
[44]. Although these models end up being very large in the number of variables
and constraints, the classical solvers available today are sophisticated and
powerful enough to handle problems of reasonable size. A significant benefit
of integer linear programming lies in its flexibility to accommodate
modifications to the problem, such as incorporating new structural motifs like
pseudoknots, without disrupting the method, unlike dynamic programming
algorithms.

Several ways have been proposed to represent the secondary structure
prediction problem, defining a variable to be:

  * •

a pair between bases $i$ and $j$ represented as $(i,j)$

  * •

a stack of pairs $(i,j),(i+1,j-1),\ldots,(i+k,j-k)$ represented as $(i,j,k)$

  * •

a stack of two consecutive pairs, also known as a quartet, represented as
$(i,j,i+1,j-1)$.

Using a pair between $i$ and $j$ as a variable works well in solving the
maximum pair problem (like Nussinov’s algorithm) but makes it difficult to
model more complex requirements and constraints. In terms of scaling, the
number of variables of the problem with respect to the sequence length scales
best when using quartets and worse when using pairs.

We use mathematical programming to model the problem as a binary optimization
problem. Each of the positions in an RNA sequence $1,\ldots,n$ consists of one
of the bases $U,A,C,G$. A pair between two bases in positions $i$ and $j$ in
the sequence is allowed if the corresponding bases form a valid pair. The
valid pairs are: $\\{(AU),(UA),(CG),(GC),(GU),(UG)\\}$. The variables are
defined for each quartet [44] which is two consecutive pairs (also referred to
as stacked pairs):

$x(i,j,i+1,j-1)=\left\\{\begin{array}[]{rl}1&\text{if $(i,j)$, $(i+1,j-1)$ are
made}\\\ 0&\text{otherwise}\end{array}\right.$

These variables are generated if the pairs are valid. The following data pre-
processing steps are performed to prepare the data for the formulation:

  * •

check if two pairs can create a quartet (can nest)

  * •

check two pairs for having the same base

  * •

checks two pairs for crossing

  * •

create a list of quartets (variables of the problem)

  * •

check two quartets for crossing

Two quartets cannot be both selected if there are crossing pairs between them.
Also, each base can be paired only with one other base. These conditions
result in constraints of the form:

$x(i,j,i+1,j-1)+x(k,\ell,k+1,\ell-1)\leq 1\,.$

Since the variables are binary, the constraint above is equivalent to
requiring that:

$x(i,j,i+1,j-1)x(k,\ell,k+1,\ell-1)=0\,,$

which can be included in the objective function as a product of the two
variables with a penalty, thus preventing both of them to be 1. We note that
depending on the algorithm used to solve the problem, we may want to keep one
(with constraints) or the other (without constraints) formulation. In our
experiments, we utilize CPLEX, which demonstrates superior performance when
explicit constraints are included in the formulation.

The objective function has the following components:

  * •

Each variable (quartet) has a free energy as given in Ref. [45].

  * •

Consecutive quartets are preferred since they provide more stability. We model
this, by adding the product of the corresponding variables in the objective
function with a reward.

  * •

Certain structures are discouraged. In this work we consider the $(UA)$ pair
ending penalty as given in Ref. [45]. We implement this by adding variables
that form such ending pair to the objective function with a penalty.

For brevity of notation, we use $q_{i}$ to denote the $i^{th}$ quartet
variable. Let $Q$ be the set of valid quartets, $QC$ the set of all pairwise
crossing quartets, $QS(q_{i})$ the set of quartets that can be stacked with
quartet $q_{i}$, and $\mathit{QUA}$ the set of stacked quartets that have a
$(UA)$ end pair. The sets $Q$, $QC$, $QS(q_{i})$ and $\mathit{QUA}$, are
determined based on the given sequence in a pre-processing step. Each quartet
$q_{i}\in Q$ is formed by two consecutive (stacked) pairs $p^{i}_{1}$ and
$p^{i}_{2}$. The MFE $e_{q_{i}}$ for quartet $q_{i}$ is the MFE of $p_{1}$
followed by $p_{2}$ [45]. Two quartets $q_{1},q_{2}\in Q$ that are stacked are
rewarded by a reward $r$. If a stack ends in a $(U,A)$ pair there is a penalty
$p$. The objective function is thus:

$\min\sum_{q_{i}\in Q}e_{q_{i}}q_{i}+r\sum_{q_{i}\in Q}{\sum_{q_{j}\in
QS(q_{i})}}q_{i}q_{j}+p\sum_{q_{i}\in
Q}\sum_{q_{j}\in\mathit{QUA}}q_{i}(1-q_{j})\,.$ (1)

The constraints of the problem implement the requirement that if two quartets
$q_{1},q_{2}\in Q$ are crossing (i.e., contain crossing pairs) they cannot be
both selected:

$q_{i}+q_{j}\leq 1\,,\quad\quad\forall q_{i},q_{j}\in QC\,.$ (2)

The resulting problem is a quadratic binary optimization problem. We use CPLEX
as a classical solver to solve this problem and create a baseline. As noted
above, for larger sequences the number of constraints in this formulation is
quite large. However, modern classical solvers, such as CPLEX, can deal with
the size of the problem. In our experiments, we generated 10,000 sequences of
length varying between 15 and 60 nucleotides, which resulted in quadratic
binary problems with up to 327 variables and up to 38,946 constraints. CPLEX
was able to find optimal solutions to all instances within a few seconds per
instance. In Fig. 3 we plot the time CPLEX takes to solve a problem, measured
in internal clock ticks, as a function of length of RNA sequences. Our dataset
contains over 600 problem instances per sequence length The box plot suggests
that the time to solve the hardest instances grows exponentially with sequence
length. As a result, in the mathematical programming formulation of the
problem, classical solvers will struggle to find a solution within a
reasonable time frame even for sequences of a few hundred nucleotides.

Figure 3: Time scaling (internal clock ticks) for CPLEX as a function of
sequence length for 10,000 random RNA sequences.

##  IV Quantum approach

###  IV-A QUBO formulation

The mathematical programming formulation described above is a quadratic
constrained binary optimization problem. The corresponding problem Hamiltonian
can be obtained by transforming the constrained problem to a QUBO problem.
Other approaches [15, 14] formulate the QUBO directly and train the model to
find the appropriate penalties and rewards for the objective function. Here we
use the Qiskit [46] available methods to convert the problem to a QUBO and
calculate the penalties for the constraints.

As previously mentioned, the constraints in Eq. 2 can be incorporated into the
objective function (Eq. 1) through the introduction of a slack variable and a
penalty term. However, due to the form of the constraints and the binary
nature of the variables, the introduction of additional slack variables is
unnecessary. After adding a penalty $t$ for the constraints, the final
objective function of the resulting QUBO takes the following form:

$\begin{array}[]{lcl}\min&&\displaystyle\sum_{q_{i}\in Q}e_{q_{i}}q_{i}\\\
&&+\displaystyle r\sum_{q_{i}\in Q}{\sum_{q_{j}\in QS(q_{i})}}q_{i}q_{j}\\\
&&+\displaystyle p\sum_{q_{i}\in
Q}\sum_{q_{j}\in\mathit{QUA}}q_{i}(1-q_{j})\\\ &&+\displaystyle
t\sum_{q_{i},q_{j}\in QC}q_{i}q_{j}\,,\\\ \text{s.t.}&&q_{i}\in\\{0,1\\}\
\forall q_{i}\in Q\,.\par\end{array}$ (3)

###  IV-B CVaR-based VQE

In this work, we use Conditional Value at Risk (CVaR)-based VQE [24, 47] to
solve the QUBO problems arising from the formulation discussed in Section
IV-A. VQE is a hybrid quantum-classical optimization algorithm which uses a
unitary $U(\boldsymbol{\theta})$, characterized by the rotation angles
$\boldsymbol{\theta}$, to prepare a parameterized quantum state
$|\psi(\boldsymbol{\theta})\rangle=U(\boldsymbol{\theta})|\boldsymbol{0}\rangle$.
The probability distribution of bitstrings, and the energies corresponding to
them, sampled after measurement of the state
$|\psi(\boldsymbol{\theta})\rangle$ is used to compute a suitable objective.
The goal is to identify the optimal parameters $\boldsymbol{\theta}$ that
minimize the objective function, typically achieved through a classical
optimization scheme. We use here a CVaR-based objective that is defined as the
average of lower $\alpha$-tail of the energy distribution of the sampled
bitstring. More formally, if the set of bitstring energies sampled and sorted
in increasing order is denoted by $\\{E_{i}\\}$ where $i\in[1..K]$, then the
CVaR($\alpha$) is defined as

$\displaystyle\frac{1}{\lceil\alpha K\rceil}\sum_{i=1}^{\lceil\alpha
K\rceil}E_{i}.$ (4)

When $\alpha=1$, CVaR($\alpha$) corresponds to the expectation value
$\langle\psi(\boldsymbol{\theta})\lvert
H_{p}\rvert\psi(\boldsymbol{\theta})\rangle$ of the problem Hamiltonian
$H_{p}$, whereas when $\alpha\to 0$, only the minimum among sampled energies
is considered. Employing an $\alpha$ value in the range $(0,1]$ is expected to
improve the convergence behavior of the classical optimizer [24, 48]. More
recently, authors in Ref. [20] indicated CVaR to be robust against quantum
hardware noise and provided bounds on the noise-free expectation values,
placed by CVaR based on the noisy samples.

###  IV-C Classical Optimizer

Recently, several studies [49, 50] have examined the performance of classical
optimizers for variational schemes, with an emphasis on robustness of such
optimizers on the noisy near-term quantum hardware. In this work, we employ
the NFT algorithm [25] as the classical optimizer which was found to be
relatively better under hardware noise [26, 51]. NFT is a gradient-free
classical optimizer that falls under the class of sequential minimal
optimization algorithms such as Rotosolve [52, 53, 54, 55, 56]. Under this
scheme, each parameter $\theta_{i}$ is optimized sequentially. NFT uses the
fact that for parameterized circuits with parameters sitting only on the
single-qubit Pauli rotations (and fixed two-qubit gates), the expectation of
an observable $\langle O\rangle$ can be expressed as an univariate sine
function of $\theta_{i}$ with three unknown coefficients. Under the scheme,
this sine function is constructed by evaluating $\langle O\rangle$ at three
different values of $\theta_{i}$, while keeping rest of the parameters fixed.
The optimal $\theta_{i}$ corresponds to the angle that minimizes the sine
function. For convergence guarantees, unbiased estimate of $\langle O\rangle$
is needed. However, numerical simulations conducted in Ref. [25] found NFT
robust under finite sampling as well. In this work, we use CVaR values as a
proxy to the expectation values. Based on the MPS-simulator VQE runs, we
notice, NFT usually achieving convergence faster when CVaR with $\alpha<1$ is
used as a proxy to the expectation ($\alpha=1$).

###  IV-D Ansatz

$|0\rangle$ & [wires=3, steps=5, style=solid, rounded corners, fill=blue!20,
inner sep=1pt, background, label style=label position=below,
yshift=-0.5cm]repeat $p$ times e^-i θ_0 Y 1 e^-i θ_3 Y
$|0\rangle$ e^-i θ_1 Y Z 1 e^-i θ_4 Y
$|0\rangle$ e^-i θ_2 Y Z e^-i θ_5 Y

Figure 4: Two-local ansatz for a three-qubit case.

We use here the “two-local” ansatz involving single-qubit Pauli-Y rotations
($e^{-i\phi Y}$) and two-qubit control-Z gates, as depicted in Fig. 4. This
configuration is compatible with the NFT optimizer. For simulations, two-qubit
gates are applied in a “pairwise” fashion composed of two layers. In the first
layer, qubit $i$ is entangled with qubit $(i+1)$, for all even values of $i$
and in the second layer qubit $i$ is entangled with qubit $(i+1)$ for all odd
values of $i$. For the hardware runs, the two-qubit gates are applied between
qubits belonging to the actual hardware sub-graphs obtained after the hardware
characterization steps described in Section V-B1.

##  V Results

###  V-A Simulation runs

Figure 5: Box plot depicting the success percentages for problem sizes between
10 and 40 qubits for CVaR-VQE runs on the MPS simulator. For each problem
size, simulations were conducted on ten distinct sequences. For each sequence,
ten independent CVaR-VQE trials were performed. The success percentage is
calculated as the number of successful runs over ten trials. The average over
different sequences is indicated by blue triangles. Figure 6: Box plot of the
optimality gap $\chi$ for problem sizes between 10 and 40 qubits for CVaR-VQE
runs on the MPS simulator. The averaging scheme is the same as described in
Fig. 5.

We perform noise-free simulation of CVaR-VQE using the matrix-product state
(MPS) [57, 58] simulator provided in Qiskit [46]. Since the computational cost
of MPS-based simulator grows exponentially with entanglement [57], we limit
our simulations to shallow-depth circuits. We use two layers ($p=2$) of the
“two-local” ansatz as shown in Fig. 4. After examining several different
values of the CVaR parameter $\alpha$ (see Section IV-B) on relatively small
problem sizes, we set it equal to 0.1 for all the simulation runs. The number
of shots $N_{\mathrm{shots}}$ used in the simulations was $2^{5}$ and $2^{10}$
for the 10- and 15-qubit problems, respectively, and $2^{13}$ for the rest of
the problems. The NFT optimizer evaluates $n_{\mathrm{eval}}=3$ sets of
circuit parameters at each iteration. Each iteration consists of
$n_{\mathrm{eval}}\times N_{\mathrm{shots}}$ circuit calls. The optimization
scheme is executed for $N_{\mathrm{iter}}$ iterations until sufficient
convergence of the CVaR values is observed, typically in the range of 100 to
200 iterations.

Because classical optimizers like NFT are prone to getting stuck in local
minima, we run $N_{\mathrm{trial}}=10$ independent trials of CVaR-VQE for each
problem instance. Each trial run starts with a set of circuit parameters
initialized randomly from a uniform distribution. We examine 10 different mRNA
sequences for each qubit size.

Fig. 5 presents a box plot illustrating the success probability across various
problem instances ranging from 10 and 40 qubits. A run is deemed successful if
the bitstring with the lowest energy, identical to the solution obtained by
CPLEX, is found within the samples collected after measuring the final quantum
state $|\psi(\boldsymbol{\theta_{f}})\rangle$, where $\boldsymbol{\theta_{f}}$
denotes the circuit parameters obtained after $N_{\mathrm{iter}}$ iterations.
The success probability $p_{\mathrm{succ}}$ is then defined as
$N_{\mathrm{succ}}/N_{\mathrm{trial}}$, where $N_{\mathrm{succ}}$ is the
number of successful trials. As evident from Fig. 5, the CVaR-VQE algorithm
demonstrates robust performance, with the average success probability in the
range from 0.4 to 1.0.

In addition to the success probability, we also use the optimality gap as a
metric to assess the quality of the solutions obtained from CVaR-VQE. The
optimality gap $\chi$ is defined as

$\displaystyle\chi=\frac{\lvert
F(\boldsymbol{\theta}_{\mathrm{f}})_{\mathrm{low}}-F_{0}\rvert}{\lvert
F_{0}\rvert}\times 100\,,$ (5)

where $F(\boldsymbol{\theta}_{\mathrm{f}})_{\mathrm{low}}$ represents the
lowest value of the objective function, as defined in Eq. 3, found among the
samples at the end of $N_{\mathrm{iter}}$ iterations. $F_{0}$ denotes the
value of the objective function of the optimal solution found by CPLEX. $\chi$
is zero when the lowest energy found by the quantum optimization matches the
one found by CPLEX, and greater than zero otherwise. The smaller the deviation
from zero, the better the quality of solutions obtained using CVaR-VQE.

Fig. 6 shows a box plot of the optimality gap $\chi_{\mathrm{avg}}$, averaged
over 10 independent trials for each sequence, for qubits sizes ranging from 10
to 40. The CVaR-VQE algorithm demonstrates high performance, achieving an
average optimality gap of less than 20$\%$.

Our results indicate that the performance of the algorithm deteriorates as the
problem size increases. For the constant-depth circuits employed in this
study, such a decline with increasing problem size is anticipated. We
hypothesize that the performance of CVaR-VQE can be further enhanced by
increasing the number of layers in the ansatz, optimizing the choice of
$\alpha$, or by incorporating other recent advancements in CVaR [59]. We
intend to investigate these aspects of the algorithm in future research.

###  V-B Hardware runs

$|0\rangle$ & H 1[2,steps=3,style=dashed,rounded corners,fill=green!20, inner
xsep=2pt,background,label style=label
position=below,anchor=north,yshift=-0.2cm]swap 1
$|0\rangle$ -1 H

Figure 7: Quantum circuit for noise characterization.

In this Section, we provide details about the hardware runs carried out on the
IBM Eagle and Heron processors, including ibm_brisbane, ibm_osaka, and
ibm_torino.

####  V-B1 Qubit Selection

(a)

(b)

Figure 8: Panel a) shows the fidelity of different couplings across the three
layers obtained during the hardware characterization scheme detailed in
Section V-B1. Panel b) shows the heavy-hex layout on ibm_brisbane and
ibm_osaka with nodes on the graph representing qubits and edges representing
qubit couplings. The couplings belonging to the three layers are color coded
as blue, red, and, green, matching the fidelity plots in a). The purple
colored nodes represent the 100 qubit chain considered for selection of qubits
on ibm_brisbane based on the fidelity values.

Due to variations in the noise characteristics across the physical chip of a
quantum computer, mapping problem qubits to the hardware sub-graphs with
smaller error rates becomes crucial. To account for the cross-talk, similar in
spirit to the layered fidelity benchmarking scheme [60], we examine fidelity
across circuit layers containing disjoint set of qubit couplings. We implement
a simple two-qubit circuit (see Fig. 7) across couplings on each layer
separately. Note, that the three layers with couplings of color blue, red, and
green in Fig. 8 cover all of the neighboring interactions present on the
quantum device. Fidelity between ideal ($00$ bitstring) and the distribution
obtained from the quantum computer is indicated in Fig. 8 for the three
layers. Next, in order to select the set of “good” qubits, we consider a
1D-chain of 100 qubits, color coded as nodes in magenta in Fig. 8 [60]. Across
the chain, we exclude qubits connected to couplings with fidelities not
meeting a set threshold. For our runs we set this fidelity threshold to 0.85.
If at the end of the process, we still need more qubits, we extend our search
to the rest of the couplings not covered by the 100-qubit chain. We employ the
described process of qubit selection on all the devices once before starting
the CVaR-VQE runs.

TABLE I: Quantum device circuit and run attributes. Max iterations column indicates the maximum number of classical iteration steps carried out by the NFT optimizer. The two-qubit Echoed Cross Resonant (ECR) gate is native to ibm_brisbane and ibm_osaka. The Control-Z (CZ) gate is native to ibm_torino. |  ibm_brisbane/ibm_osaka | ibm_torino
---|---|---
Qubits |  26 |  40 |  50 | 80
Circuit depth |  18 |  20 |  20 | 17
Max iterations |  104 |  320 |  600 | –
ECR$/$CZ Count |  25 |  39 |  49 | 158

(a) 26 qubits, 25 nucleotides

(b) 40 qubits, 30 nucleotides

(c) 50 qubits, 30 nucleotides

Figure 9: Figures show actual optimal sequences based on the lowest energy
bitstring found during hardware runs.

####  V-B2 Error suppression and error mitigation

We use Dynamical Decoupling (DD) [61, 62, 63, 64], which adds pulses during
the idle times of the circuit and helps to suppress errors without any
additional overhead. Additionally, we use a version of Matrix-free Measurement
Mitigation (M3) [65] for readout error mitigation. Instead of working based
off a full assignment matrix of size $2^{n}$, M3 uses a much smaller subspace
formed by the noisy bitstrings, allowing it to usually handle large system
efficiently [65]. In order to keep the error mitigation overhead to “minimal”,
we do not consider mitigation schemes like ZNE [23, 66, 67] and PEC [68] in
this work. Interestingly, as pointed out in Ref. [20], CVaR-based schemes with
a relatively small sampling overhead may still be able to track the noise-free
expectation values [20]. We plan to explore these aspects of the CVaR
algorithm in near future. Additionally, we set optimization_level=3 under the
Qiskit’s Sampler V1 primitive runtime options to optimize circuits, which
gives rise to transpiled circuits depths in the range of 17-20 (see Table I).

Figure 10: Success probability for different number of qubits for the hardware
runs. Figure 11: Box plot of the optimality gap $\chi$ for problem sizes in
the range of 10 to 40 qubits executed on quantum hardware. Average $\chi$ over
$N_{\mathrm{trial}}$ independent runs is shown using blue triangles.

####  V-B3 Hardware Results

On the hardware, we study 3 mRNA sequences: (a) of length 25, represented by
the problem instance with 26 qubits, (b) of length 30, with 40 qubits, and (c)
of length 30, with 50 qubits. We run 8, 8, and 5 independent CVaR-VQE trials
for these instances, correspondingly, using $\alpha=0.25$ and $p=1$.

The success probability of CVaR-VQE experiments on quantum hardware is shown
in Fig. 10. It should be noted that hardware experiments were conducted for
only one sequence of each problem size; therefore, the success probability is
computed as a cumulative quantity over 10 trials. As indicated by the plot,
the success probability across different qubit sizes is above $\sim 0.38$. The
average optimality gap $\chi_{\mathrm{avg}}$ of CVaR-VQE experiments on the
hardware is shown in Fig. 11 and is below $\sim$34$\%$. We set the number of
shots for each circuit to $2^{8}$ for the 26-qubit problem and $2^{13}$ for
the larger problem instances. The maximum number of iterations required to
find the optimal parameters on the hardware is reported in the Table I. Some
trials took fewer than the maximum number of iterations reported. The lowest
energy sequences found are shown in Fig. 9.

Figure 12: (a) Measurement probability of sampled bitstrings plotted against
bitstring energies for the 80-qubit problem. Red bar indicates the probability
corresponding to the lowest energy bitstring. (b) Optimal folded structure of
the 42-nucleotide, 80-qubit mRNA sequence based on the corresponding lowest
energy bitstring found by the hardware run.

For the 42-nucleotide, 80-qubit mRNA problem instance, we do not carry out
circuit parameter optimization on the hardware. The optimal parameters for the
circuit were obtained by running CVaR-VQE using the MPS simulator with the
$\alpha=0.1$ and $p=2$. We use the optimized parameters to run circuits on the
recently introduced ibm_torino device featuring the 133-qubit tunable-coupler
Heron processor [69]. The resulting circuit has depth 17 and contains 158 two-
qubit CZ gates. We plot the distribution of measurement probabilities of the
sampled bitstrings with the corresponding energies in Fig. 12. The plot
illustrates that the lowest energy bitstring with energy value of -140, which
matches the lowest energy bitstring found by CPLEX, is sampled with a
probability separated from the rest of the sampling probabilities by more than
one order of magnitude. The measurement probability of the lowest bitstring is
$0.045$, hence, we need at most $1/0.045\sim 20$ samples to find the lowest
energy bitstring. The optimal sequence found is shown in Fig. 12. Note, the
offline training of circuit parameters using MPS simulator was only possible
because of the classical simulatability of such circuits. As pointed out in
Refs. [21, 70, 71, 72, 73, 74, 75, 76, 77], the heavy-hex lattice-based
circuits characterized by depths we employ in this study can be efficiently
simulated using classical schemes. Our results point to an exciting regime
where relatively large scale optimization-based problems, approaching utility-
scale, can also be successfully executed on the noisy quantum hardware and can
produce reliable results using low-overhead error mitigation schemes.

###  V-C Feasibility of hardware runs

Unlike schemes targeting fault-tolerant quantum devices, where quantum
resource scaling arguments are relatively straightforward to derive,
establishing such arguments for variational schemes is more challenging [78].
If $p$ layers of an ansatz are employed, and each circuit layer requires time
$T$ to execute, the total time taken is approximately
$N_{\mathrm{circuit}}\times p\times T$, where $N_{\mathrm{circuit}}$
represents the total number of circuits executed. Given that this study
utilizes fixed-depth circuits and assuming the execution time $T$ remains
relatively invariant with problem size, the primary factor that necessitates
closer examination is $N_{\mathrm{circuit}}$.

The number of circuits $N_{\mathrm{circuit}}$ executed during the scheme is
determined by the expression $1/p_{\mathrm{succ}}\times
N_{\mathrm{iter}}\times N_{\mathrm{shots}}\times n_{\mathrm{eval}}$, where
$n_{\mathrm{eval}}$ represents the number of different sets of circuit
parameters evaluated at each iteration to suggest a new set of circuit
parameters $\boldsymbol{\theta}$111For the NFT optimizer,
$n_{\mathrm{eval}}=3$. To assess the feasibility of solving structure
prediction problems on a quantum computer, it is crucial that none of the
three factors
$\\{1/p_{\mathrm{succ}},N_{\mathrm{iter}},N_{\mathrm{shots}}\\}$, with
$n_{\mathrm{eval}}$ assumed constant, scale exponentially with the number of
qubits. A concerning trend observed is the decline in success probabilities
$p_{\mathrm{succ}}$ as the problem size increases. However, based on the
results, it is difficult to definitively conclude that the decline is
exponential. Additionally, since the choice of algorithm dictates the lower
bound on $p_{\mathrm{succ}}$, there is potential for improvement through
algorithmic advances. Unfortunately, beyond a certain regime of classical
simulatability, performance comparison of different optimization schemes can
only be carried out in conjunction with the quantum hardware.

##  VI Conclusions and outlook

In this work, we examined the feasibility of solving secondary structure
prediction-based optimization problems of RNA sequences on a quantum computer.
Our results demonstrate that optimization schemes like CVaR-VQE running on
noisy quantum hardware with careful calibration and nominal error-mitigation,
can produce sufficiently reliable outputs even for relatively large problem
instances. However, the scalability of these methods to larger problem sizes
exceeding 100 qubits remains an open question that warrants further
investigation. As quantum hardware continues to advance, any improvements in
algorithmic techniques will further bolster prospects of using quantum
computers for optimization. This study, which primarily focused on the
capabilities of utility-scale quantum hardware, only considered a subset of
available quantum approaches; however, our future work will explore more
advanced optimization techniques.

Our results indicate that quantum computing offers valuable insights into RNA
folding mechanisms and enables rapid identification of functional RNA
structures. This work lays the groundwork for future research in quantum-
assisted bioinformatics and highlights the importance of developing scalable
quantum algorithms for real-world biological problems.

## Acknowledgments

The authors would like to thank Yashrajsinh Jadeja and Haining Lin for their
valuable assistance in acquiring test sequences and for their contributions to
the prediction of classical RNA secondary structures.

## References

  * [1] Q. Vicens and J. S. Kieft, “Thoughts on how to think (and talk) about RNA structure,” _Proceedings of the National Academy of Sciences_ , vol. 119, no. 17, p. e2112677119, 2022.
  * [2] M. Metkar, C. S. Pepin, and M. J. Moore, “Tailor made: the art of therapeutic mRNA design,” _Nature Reviews Drug Discovery_ , vol. 23, no. 1, pp. 67–83, 2024.
  * [3] S. Qin, X. Tang, Y. Chen, K. Chen, N. Fan, W. Xiao, Q. Zheng, G. Li, Y. Teng, M. Wu _et al._ , “mRNA-based therapeutics: powerful and versatile tools to combat diseases,” _Signal Transduction and Targeted Therapy_ , vol. 7, no. 1, p. 166, 2022.
  * [4] D. H. Mathews, W. N. Moss, and D. H. Turner, “Folding and finding RNA secondary structure,” _Cold Spring Harbor Perspectives in Biology_ , vol. 2, no. 12, p. a003665, 2010.
  * [5] R. B. Lyngsø and C. N. Pedersen, “RNA pseudoknot prediction in energy-based models,” _Journal of Computational Biology_ , vol. 7, no. 3-4, pp. 409–427, 2000.
  * [6] M. Zuker, “Computer prediction of RNA structure,” in _Methods in enzymology_. Elsevier, 1989, vol. 180, pp. 262–288.
  * [7] M. Zuker, “Mfold web server for nucleic acid folding and hybridization prediction,” _Nucleic Acids Research_ , vol. 31, no. 13, pp. 3406–3415, 2003.
  * [8] P. Schuster, P. F. Stadler, and A. Renner, “RNA structures and folding: from conventional to new issues in structure predictions,” _Current Opinion in Structural Biology_ , vol. 7, no. 2, pp. 229–235, 1997.
  * [9] M. Zuker, “On finding all suboptimal foldings of an RNA molecule,” _Science_ , vol. 244, no. 4900, pp. 48–52, 1989.
  * [10] E. P. Nawrocki and S. R. Eddy, “Infernal 1.1: 100-fold faster RNA homology searches,” _Bioinformatics_ , vol. 29, no. 22, pp. 2933–2935, 2013.
  * [11] S. Zakov, Y. Goldberg, M. Elhadad, and M. Ziv-Ukelson, “Rich parameterization improves RNA structure prediction,” _Journal of Computational Biology_ , vol. 18, no. 11, pp. 1525–1542, 2011.
  * [12] M. Szikszai, M. Wise, A. Datta, M. Ward, and D. H. Mathews, “Deep learning models for RNA secondary structure prediction (probably) do not generalize across families,” _Bioinformatics_ , vol. 38, no. 16, pp. 3892–3899, 2022.
  * [13] H. Zhang, L. Zhang, D. H. Mathews, and L. Huang, “Linearpartition: linear-time approximation of RNA folding partition function and base-pairing probabilities,” _Bioinformatics_ , vol. 36, no. Supplement_1, pp. i258–i267, 2020.
  * [14] D. M. Fox, C. M. MacDermaid, A. M. Schreij, M. Zwierzyna, and R. C. Walker, “RNA folding using quantum computers,” _PLOS Computational Biology_ , vol. 18, no. 4, p. e1010032, 2022.
  * [15] T. Zaborniak, J. Giraldo, H. Müller, H. Jabbari, and U. Stege, “A QUBO model of the RNA folding problem optimized by variational hybrid quantum annealing,” in _2022 IEEE International Conference on Quantum Computing and Engineering (QCE)_. IEEE, 2022, pp. 174–185.
  * [16] J. Jiang, Q. Yan, Y. Li, Y. Chai, M. Lu, Z. Cui, M. Dou, Q. Wang, Y.-C. Wu, and G.-P. Guo, “Predicting RNA secondary structure on universal quantum computer,” _arXiv preprint arXiv:2305.09561_ , 2023.
  * [17] E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimization algorithm,” _arXiv preprint arXiv:1411.4028_ , 2014.
  * [18] M. P. Harrigan, K. J. Sung, M. Neeley, K. J. Satzinger, F. Arute, K. Arya, J. Atalaya, J. C. Bardin, R. Barends, S. Boixo _et al._ , “Quantum approximate optimization of non-planar graph problems on a planar superconducting processor,” _Nature Physics_ , vol. 17, no. 3, pp. 332–336, 2021.
  * [19] Z. He, R. Shaydulin, S. Chakrabarti, D. Herman, C. Li, Y. Sun, and M. Pistoia, “Alignment between initial state and mixer improves QAOA performance for constrained optimization,” _npj Quantum Information_ , vol. 9, no. 1, p. 121, 2023.
  * [20] S. V. Barron, D. J. Egger, E. Pelofske, A. Bärtschi, S. Eidenbenz, M. Lehmkuehler, and S. Woerner, “Provable bounds for noise-free expectation values computed from noisy samples,” _arXiv preprint arXiv:2312.00733_ , 2023.
  * [21] S. H. Sack and D. J. Egger, “Large-scale quantum approximate optimization on nonplanar graphs with machine learning noise mitigation,” _Physical Review Research_ , vol. 6, no. 1, p. 013223, 2024.
  * [22] A. Miessen, D. J. Egger, I. Tavernelli, and G. Mazzola, “Benchmarking digital quantum simulations and optimization above hundreds of qubits using quantum critical dynamics,” _arXiv preprint arXiv:2404.08053_ , 2024.
  * [23] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. Van Den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme _et al._ , “Evidence for the utility of quantum computing before fault tolerance,” _Nature_ , vol. 618, no. 7965, pp. 500–505, 2023.
  * [24] P. K. Barkoutsos, G. Nannicini, A. Robert, I. Tavernelli, and S. Woerner, “Improving Variational Quantum Optimization using CVaR,” _Quantum_ , vol. 4, p. 256, Apr. 2020. doi: 10.22331/q-2020-04-20-256. [Online]. Available: https://doi.org/10.22331/q-2020-04-20-256
  * [25] K. M. Nakanishi, K. Fujii, and S. Todo, “Sequential minimal optimization for quantum-classical hybrid algorithms,” _Phys. Rev. Res._ , vol. 2, p. 043158, Oct 2020. doi: 10.1103/PhysRevResearch.2.043158
  * [26] M. Oliv, A. Matic, T. Messerer, and J. M. Lorenz, “Evaluating the impact of noise on the performance of the variational quantum eigensolver,” _arXiv preprint arXiv:2209.12803_ , 2022.
  * [27] IBM, “IBM ILOG CPLEX optimizer
https://www.ibm.com/products/ilog-cplex-optimization-studio/cplex-optimizer.”

  * [28] A. Mamuye, E. Merelli, and L. Tesei, “A graph grammar for modelling RNA folding,” _Electronic Proceedings in Theoretical Computer Science_ , vol. 231, pp. 31–41, 12 2016. doi: 10.4204/EPTCS.231.3
  * [29] ViennaRNA, “ViennaRNA server http://rna.tbi.univie.ac.at/.”
  * [30] R. Durbin, S. R. Eddy, A. Krogh, and G. Mitchison, _Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids_. Cambridge university press, 1998.
  * [31] M. Andronescu, A. Condon, H. H. Hoos, D. H. Mathews, and K. P. Murphy, “Efficient parameter estimation for RNA secondary structure prediction,” _Bioinformatics_ , vol. 23, no. 13, pp. i19–i28, 2007.
  * [32] J. Zuber, S. J. Schroeder, H. Sun, D. H. Turner, and D. H. Mathews, “Nearest neighbor rules for RNA helix folding thermodynamics: improved end effects,” _Nucleic Acids Research_ , vol. 50, no. 9, pp. 5251–5262, 2022.
  * [33] M. W. Lewis, A. Verma, and T. T. Eckdahl, “Qfold: a new modeling paradigm for the RNA folding problem,” _Journal of Heuristics_ , vol. 27, pp. 695–717, 2021. [Online]. Available: https://doi.org/10.1007/s10732-021-09471-3
  * [34] J.-H. Chen, S.-Y. Le, and J. V. Maizel, “Prediction of common secondary structures of RNAs: a genetic algorithm approach,” _Nucleic Acids Research_ , vol. 28, no. 4, pp. 991–999, 02 2000. doi: 10.1093/nar/28.4.991. [Online]. Available: https://doi.org/10.1093/nar/28.4.991
  * [35] Q. Zhao, Z. Zhao, X. Fan, Z. Yuan, Q. Mao, and Y. Yao, “Review of machine learning methods for RNA secondary structure prediction,” _PLoS Comput Biol_ , vol. 17, no. 8: e1009291, 2021. [Online]. Available: https://doi.org/10.1371/journal.pcbi.1009291
  * [36] D. H. Seetin, Matthew G.and Mathews, _RNA Structure Prediction: An Overview of Methods_. Totowa, NJ: Humana Press, 2012, pp. 99–122. ISBN 978-1-61779-949-5. [Online]. Available: https://doi.org/10.1007/978-1-61779-949-5_8
  * [37] S. R. Eddy, “How do RNA folding algorithms work?” _Nature Biotechnology_ , pp. 1457–1458, 2004. doi: 10.1038/nbt1104-1457. [Online]. Available: https://doi.org/10.1038/nbt1104-1457
  * [38] M. Zuker and P. Stiegler, “Optimal computer folding of large RNA sequences using thermodynamics and auxiliary information,” _Nucleic Acids Research_ , vol. 9, no. 1, pp. 133–148, 1981.
  * [39] D. Draper, T. Gluick, and P. Schlax, “In RNA structure and function, vol. 298,” 2000.
  * [40] S. D. Gilbert, R. P. Rambo, D. Van Tyne, and R. T. Batey, “Structure of the SAM-II riboswitch bound to s-adenosylmethionine,” _Nature Structural & Molecular Biology_, vol. 15, no. 2, pp. 177–182, 2008.
  * [41] D. J. Klein, T. E. Edwards, and A. R. Ferré-D’Amaré, “Cocrystal structure of a class I preQ1 riboswitch reveals a pseudoknot recognizing an essential hypermodified nucleobase,” _Nature Structural & Molecular Biology_, vol. 16, no. 3, pp. 343–344, 2009.
  * [42] R. C. Spitale, A. T. Torelli, J. Krucinska, V. Bandarian, and J. E. Wedekind, “The structural basis for recognition of the PreQ0 metabolite by an unusually small riboswitch aptamer domain*,” _Journal of Biological Chemistry_ , vol. 284, no. 17, pp. 11 012–11 016, 2009.
  * [43] D. H. Mathews, J. Sabina, M. Zuker, and D. H. Turner, “Expanded sequence dependence of thermodynamic parameters improves prediction of RNA secondary structure,” _Journal of Molecular Biology_ , vol. 6, pp. 911–940, 1999.
  * [44] D. Gusfield, _Integer Linear Programming in Computational and Systems Biology: An Entry-Level Text and Course_. Cambridge University Press, 2019.
  * [45] D. H. Turner and D. H. Mathews, “NNDB: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure,” _Nucleic Acids Research_ , vol. 38, no. suppl_1, p. D280– D282, 2009.
  * [46] Qiskit contributors, “Qiskit: An open-source framework for quantum computing,” 2023.
  * [47] A. Robert, P. K. Barkoutsos, S. Woerner, and I. Tavernelli, “Resource-efficient quantum algorithm for protein folding,” _npj Quantum Information_ , vol. 7, no. 1, p. 38, 2021.
  * [48] I. Kolotouros and P. Wallden, “An evolving objective function for improved variational quantum optimisation,” _arXiv preprint arXiv:2105.11766_ , 2021.
  * [49] W. Lavrijsen, A. Tudor, J. Müller, C. Iancu, and W. De Jong, “Classical optimizers for noisy intermediate-scale quantum devices,” in _2020 IEEE International Conference on Quantum Computing and Engineering (QCE)_. IEEE, 2020, pp. 267–277.
  * [50] A. Pellow-Jarman, S. McFarthing, I. Sinayskiy, A. Pillay, and F. Petruccione, “QAOA performance in noisy devices: the effect of classical optimizers and ansatz depth,” _arXiv preprint arXiv:2307.10149_ , 2023.
  * [51] H. Singh, S. Majumder, and S. Mishra, “Benchmarking of different optimizers in the variational quantum algorithms for applications in quantum chemistry,” _The Journal of Chemical Physics_ , vol. 159, no. 4, 2023.
  * [52] M. Ostaszewski, E. Grant, and M. Benedetti, “Structure optimization for parameterized quantum circuits,” _Quantum_ , vol. 5, p. 391, 2021.
  * [53] W. Li, Y. Ge, S.-X. Zhang, Y.-Q. Chen, and S. Zhang, “Efficient and robust parameter optimization of the unitary coupled-cluster ansatz,” _Journal of Chemical Theory and Computation_ , 2024.
  * [54] J. G. Vidal and D. O. Theis, “Calculus on parameterized quantum circuits,” _arXiv preprint arXiv:1812.06323_ , 2018.
  * [55] R. M. Parrish, J. T. Iosue, A. Ozaeta, and P. L. McMahon, “A jacobi diagonalization and anderson acceleration algorithm for variational quantum algorithm parameter optimization,” _arXiv preprint arXiv:1904.03206_ , 2019.
  * [56] D. Wierichs, J. Izaac, C. Wang, and C. Y.-Y. Lin, “General parameter-shift rules for quantum gradients,” _Quantum_ , vol. 6, p. 677, 2022.
  * [57] G. Vidal, “Efficient classical simulation of slightly entangled quantum computations,” _Physical Review Letters_ , vol. 91, no. 14, p. 147902, 2003.
  * [58] U. Schollwöck, “The density-matrix renormalization group in the age of matrix product states,” _Annals of Physics_ , vol. 326, no. 1, pp. 96–192, 2011.
  * [59] I. Kolotouros and P. Wallden, “Evolving objective function for improved variational quantum optimization,” _Physical Review Research_ , vol. 4, no. 2, p. 023225, 2022.
  * [60] D. C. McKay, I. Hincks, E. J. Pritchett, M. Carroll, L. C. Govia, and S. T. Merkel, “Benchmarking quantum processor performance at scale,” _arXiv preprint arXiv:2311.05933_ , 2023.
  * [61] L. Viola and S. Lloyd, “Dynamical suppression of decoherence in two-state quantum systems,” _Physical Review A_ , vol. 58, no. 4, p. 2733, 1998.
  * [62] P. Zanardi, “Symmetrizing evolutions,” _Physics Letters A_ , vol. 258, no. 2-3, pp. 77–82, 1999.
  * [63] D. Vitali and P. Tombesi, “Using parity kicks for decoherence control,” _Physical Review A_ , vol. 59, no. 6, p. 4178, 1999.
  * [64] N. Ezzell, B. Pokharel, L. Tewala, G. Quiroz, and D. A. Lidar, “Dynamical decoupling for superconducting qubits: a performance survey,” _Physical Review Applied_ , vol. 20, no. 6, p. 064027, 2023.
  * [65] P. D. Nation, H. Kang, N. Sundaresan, and J. M. Gambetta, “Scalable mitigation of measurement errors on quantum computers,” _PRX Quantum_ , vol. 2, no. 4, p. 040326, 2021.
  * [66] A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, “Error mitigation extends the computational reach of a noisy quantum processor,” _Nature_ , vol. 567, no. 7749, pp. 491–495, 2019.
  * [67] Y. Li and S. C. Benjamin, “Efficient variational quantum simulator incorporating active error minimization,” _Physical Review X_ , vol. 7, no. 2, p. 021050, 2017.
  * [68] E. Van Den Berg, Z. K. Minev, A. Kandala, and K. Temme, “Probabilistic error cancellation with sparse Pauli–Lindblad models on noisy quantum processors,” _Nature Physics_ , vol. 19, no. 8, pp. 1116–1121, 2023.
  * [69] IBM Newsroom. IBM debuts next-generation quantum processor & IBM quantum system two, extends roadmap to advance era of quantum utility. [Online]. Available: https://newsroom.ibm.com/2023-12-04-IBM-Debuts-Next-Generation-Quantum-Processor-IBM-Quantum-System-Two,-Extends-Roadmap-to-Advance-Era-of-Quantum-Utility
  * [70] S. Patra, S. S. Jahromi, S. Singh, and R. Orús, “Efficient tensor network simulation of IBM’s largest quantum processors,” _Physical Review Research_ , vol. 6, no. 1, p. 013326, 2024.
  * [71] J. Tindall, M. Fishman, E. M. Stoudenmire, and D. Sels, “Efficient tensor network simulation of IBM’s Eagle kicked Ising experiment,” _PRX Quantum_ , vol. 5, no. 1, p. 010308, 2024.
  * [72] H.-J. Liao, K. Wang, Z.-S. Zhou, P. Zhang, and T. Xiang, “Simulation of IBM’s kicked Ising experiment with projected entangled pair operator,” _arXiv preprint arXiv:2308.03082_ , 2023.
  * [73] S. Anand, K. Temme, A. Kandala, and M. Zaletel, “Classical benchmarking of zero noise extrapolation beyond the exactly-verifiable regime,” _arXiv preprint arXiv:2306.17839_ , 2023.
  * [74] T. Begušić and G. K. Chan, “Fast classical simulation of evidence for the utility of quantum computing before fault tolerance,” _arXiv preprint arXiv:2306.16372_ , 2023.
  * [75] T. Begušić, J. Gray, and G. K.-L. Chan, “Fast and converged classical simulations of evidence for the utility of quantum computing before fault tolerance,” _Science Advances_ , vol. 10, no. 3, p. eadk4321, 2024.
  * [76] M. S. Rudolph, E. Fontana, Z. Holmes, and L. Cincio, “Classical surrogate simulation of quantum systems with lowesa,” _arXiv preprint arXiv:2308.09109_ , 2023.
  * [77] Y. Shao, F. Wei, S. Cheng, and Z. Liu, “Simulating quantum mean values in noisy variational quantum algorithms: A polynomial-scale approach,” _arXiv preprint arXiv:2306.05804_ , 2023.
  * [78] G. Scriva, N. Astrakhantsev, S. Pilati, and G. Mazzola, “Challenges of variational quantum optimization with measurement shot noise,” _Physical Review A_ , vol. 109, no. 3, p. 032408, 2024.