Resource-Efficient Quantum Algorithm for Protein Folding

TL;DR

This paper introduces a resource-efficient quantum algorithm for protein folding, using a variational approach on NISQ devices to simulate small peptide structures with promising accuracy and scalability.

Contribution

It presents a novel quantum variational algorithm with an $ ext{O}(N^4)$ model Hamiltonian for protein folding, bridging coarse-grained models and lattice representations on NISQ hardware.

Findings

Successfully simulated 10 amino acid peptide folding on 22 qubits.

Applied the method to a 7 amino acid neuropeptide on a 9-qubit IBM Q device.

Demonstrated the potential of hybrid quantum-classical algorithms for biological problems.

Abstract

Predicting the three-dimensional (3D) structure of a protein from its primary sequence of amino acids is known as the protein folding (PF) problem. Due to the central role of proteins' 3D structures in chemistry, biology and medicine applications (e.g., in drug discovery) this subject has been intensively studied for over half a century. Although classical algorithms provide practical solutions, sampling the conformation space of small proteins, they cannot tackle the intrinsic NP-hard complexity of the problem, even reduced to its simplest Hydrophobic-Polar model. While fault-tolerant quantum computers are still beyond reach for state-of-the-art quantum technologies, there is evidence that quantum algorithms can be successfully used on Noisy Intermediate-Scale Quantum (NISQ) computers to accelerate energy optimization in frustrated systems. In this work, we present a model Hamiltonian with $\mathcal{O}(N^4)$ scaling and a corresponding quantum variational algorithm for the folding of a polymer chain with $N$ monomers on a tetrahedral lattice. The model reflects many physico-chemical properties of the protein, reducing the gap between coarse-grained representations and mere lattice models. We use a robust and versatile optimisation scheme, bringing together variational quantum algorithms specifically adapted to classical cost functions and evolutionary strategies (genetic algorithms), to simulate the folding of the 10 amino acid Angiotensin peptide on 22 qubits. The same method is also successfully applied to the study of the folding of a 7 amino acid neuropeptide using 9 qubits on an IBM Q 20-qubit quantum computer. Bringing together recent advances in building gate-based quantum computers with noise-tolerant hybrid quantum-classical algorithms, this work paves the way towards accessible and relevant scientific experiments on real quantum processors.

Tables

Table 1. Table S1: Value of the turn parameter t i ( 1 ) subscript 𝑡 superscript 𝑖 1 t_{i^{(1)}} imposed by chirality constraint as a function of the values t i − 1 subscript 𝑡 𝑖 1 t_{i-1} and t i subscript 𝑡 𝑖 t_{i} for the case i 𝑖 i even (i.e., belonging the the sub-lattice 𝒜 𝒜 \mathcal{A} ). The case i 𝑖 i odd is identified by transposing the matrix.

	0	$1$	$2$	$3$
$\overline{0}$		$2$	$3$	$1$
$\overline{1}$	$3$		$0$	$2$
$\overline{2}$	$1$	$3$		$0$
$\overline{3}$	$2$	$0$	$1$

Table 2. Table S2: Comparison of existing models for 3D protein native structure prediction using quantum algorithms. The sore point of the existing models (Perdomo  et al. , 2008 ; Babbush  et al. , 2012 ; Perdomo-Ortiz  et al. , 2012 ; Babej  et al. , 2018 ; Fingerhuth  et al. , 2018 ) is the locality of the generated Hamiltonian and/or the required number of ancillas (qubit resources). In ‘digital’ quantum computers (last two columns), the coupling between more than two qubits can be decomposed into two-qubit controlled operations at the cost of an increased circuit depth. On the other hand, the architecture of analog quantum annealers prevents the optimization of non-Ising like Hamiltonian (i.e., with k 𝑘 k -body interactions between qubits). In fact, locality-reduction to an isospectral Ising Hamiltonian ( k = 2 𝑘 2 k=2 ) requires the introduction of additional ancillas qubits and the use of complex mapping schemes  (Kempe  et al. , 2006 ; Barkoutsos  et al. , 2018 ) (first three columns). Since k 𝑘 k -local terms require k − 2 𝑘 2 k-2 ancillas to be embedded, the resource requirements for both quantum approaches (digital and analog) scales exponentially  (Babbush  et al. , 2012 ) with the polymer size if the locality of the model depends on the polymer length, N 𝑁 N , for coarse-grained models. While for the model presented in this work the locality is independent from N 𝑁 N this is not the case for most of the proposed Hamiltonians in the literature (Perdomo  et al.   (Perdomo  et al. , 2008 ) , Perdomo  et al.   (Perdomo-Ortiz  et al. , 2012 ) , and Babej  et al.   (Babej  et al. , 2018 ; Fingerhuth  et al. , 2018 ) ). The corresponding scalings and the localities are defined in the text. For our model, they are obtained after constructing the Hamiltonian on Qiskit  (Aleksandrowicz  et al. , 2019 ) . The values recorded are obtained without use of any locality reduction schemes.

	Perdomo et al. (Perdomo et al., 2008)	Perdomo et al. (Perdomo-Ortiz et al., 2012)	Babbush et al. (Babbush et al., 2012)	Babej et al. (Babej et al., 2018; Fingerhuth et al., 2018)	This model
Model	Hydrophobic-Polar (HP)	Coarse-Grained	Coarse-Grained	Coarse-Grained	Coarse-Grained
Lattice	all	all	all	all	Tetrahedral
Types	$2$	$\infty$	$\infty$	$\infty$	$\infty$
Interactions	Nearest	Nearest	Nearest	Nearest	$l^{\mathrm{th}}$ Nearest
Locality	${\log }_{2}N$	$N$	$4$	$N$	$l+2$
Qubits	$N\log \left(N\right)$	$N^2\log \left(N\right)$	$N^3\log \left(N\right)$	$N$	$N^2\exp \left(l\right)$
Scaling	$N^8$	$\exp \left(N\right)$	$N^{12}{\log }^4\left(N\right)$	$\exp \left(N\right)$	$N^4$
Experiment	No	D-Wave	No	Rigetti QPU	IBM QPU

Equations

$$H(\mathrm{q})=H_{\text{gc}}({\mathrm{q}}_{\text{cf}})+H_{\text{ch}}({\mathrm{q}}_{\text{cf}})+H_{\text{in}}(\mathrm{q}).$$

$$H(\mathrm{q})=\sum _{\mathbf{\gamma }}^{N_{\mathrm{t}}}{h}_{\mathbf{\gamma }}\underset{i=1}{\overset{n}{⨂}}{\mathrm{q}}_i^{{\gamma }_i}$$

$$\left[0010\right]\left[0001({\mathrm{q}}_{2^{(1)}}^{(1)}{\mathrm{q}}_{2^{(1)}}^{(2)}{\mathrm{q}}_{2^{(1)}}^{(3)}{\mathrm{q}}_{2^{(1)}}^{(4)}...)\right]\left[{\mathrm{q}}_9{\mathrm{q}}_{10}{\mathrm{q}}_{11}{\mathrm{q}}_{12}({\mathrm{q}}_{3^{(1)}}^{(1)}{\mathrm{q}}_{3^{(1)}}^{(2)}{\mathrm{q}}_{3^{(1)}}^{(3)}{\mathrm{q}}_{3^{(1)}}^{(4)}...)\right]...\left[{\mathrm{q}}_{4(N-1)-3}{\mathrm{q}}_{4(N-1)-2}{\mathrm{q}}_{4(N-1)-1}{\mathrm{q}}_{4(N-1)}\right]$$

$$\left[01\right]\left[00({\mathrm{q}}_{2^{(1)}}^{(1)}{\mathrm{q}}_{2^{(1)}}^{(2)}...)\right]\left[{\mathrm{q}}_5{\mathrm{q}}_6({\mathrm{q}}_{3^{(1)}}^{(1)}{\mathrm{q}}_{3^{(1)}}^{(2)}...)\right]\left[{\mathrm{q}}_7{\mathrm{q}}_8({\mathrm{q}}_{4^{(1)}}^{(1)}{\mathrm{q}}_{4^{(1)}}^{(2)}...)\right]...\left[{\mathrm{q}}_{2(N-1)-1}{\mathrm{q}}_{2(N-1)}\right]$$

$$f_0(i)=(1-{\mathrm{q}}_{2i-1})(1-{\mathrm{q}}_{2i})$$

$$f_1(i)={\mathrm{q}}_{2i}({\mathrm{q}}_{2i}-{\mathrm{q}}_{2i-1})$$

$$f_2(i)={\mathrm{q}}_{2i-1}({\mathrm{q}}_{2i-1}-{\mathrm{q}}_{2i})$$

$$f_3(i)={\mathrm{q}}_{2i-1}{\mathrm{q}}_{2i}.$$

$$f_0(i)={\mathrm{q}}_{4i-3}$$

$$f_1(i)={\mathrm{q}}_{4i-2}$$

$$f_2(i)={\mathrm{q}}_{4i-1}$$

$$f_3(i)={\mathrm{q}}_{4i}$$

$$\Delta n_a(i,j)=\sum _{k=i}^{j-1}(-1{)}^k{f}_a(k)$$

$$\Delta n_a(i^{(s)},j^{(p)})=\Delta n_a(i,j)+\sum _{l=1}^p(-1{)}^l{f}_a(j^{(l)})-\sum _{m=1}^s(-1{)}^m{f}_a(i^{(m)})$$

$$\mathbf{x}(i,j)=\left(\begin{array}{c}\Delta n_1(i,j)\\ \Delta n_2(i,j)\\ \Delta n_3(i,j)\\ \Delta n_4(i,j)\end{array}\right)$$

$$\begin{array}[]{cccc}d(i,j)=0\Rightarrow r_{ij}=0\\ d(i,j)=1\Rightarrow r_{ij}=1\\ d(i,j)=2\Rightarrow r_{ij}=2\sqrt{\frac{2}{3}}\simeq 1.63\\ d(i,j)=3\Rightarrow r_{ij}=\sqrt{\frac{11}{3}}\simeq 1.91\\ d(i,j)=4\Rightarrow r_{ij}=\frac{4}{\sqrt{3}}\simeq 2.31\\ d(i,j)=5\Rightarrow r_{ij}=\sqrt{\frac{19}{3}}\simeq 2.52\\ ...&\end{array}$$

$$H_{\text{optional}}=\sum _{i=3}^{N-1}\lambda ({\mathrm{q}}_{4i-3}+{\mathrm{q}}_{4i-2}+{\mathrm{q}}_{4i-1}+{\mathrm{q}}_{4i}-1{)}^2$$

$$T(i,j)=\sum _{a=\{0,1,2,3\}}{f}_a(i)f_a(j)$$

$$H_{\mathrm{gc}}=\sum _{i=3}^{N-1}{\lambda }_{\text{back}}T(i,i+1)$$

$$f_0^{\text{ex}}(i^{(1)})$$

$$\displaystyle+g_{i}\Bigl{(}f_{1}(i)f_{2}(i-1)+f_{2}(i)f_{3}(i-1)+f_{3}(i)f_{1}(i-1)\Bigr{)}$$

$$f_1^{\text{ex}}(i^{(1)})$$

$$\displaystyle+g_{i}\Bigl{(}f_{0}(i)f_{3}(i-1)+f_{2}(i)f_{0}(i-1)+f_{3}(i)f_{2}(i-1)\Bigr{)}$$

$$f_2^{\text{ex}}(i^{(1)})$$

$$\displaystyle+g_{i}\Bigl{(}f_{0}(i)f_{1}(i-1)+f_{1}(i)f_{3}(i-1)+f_{3}(i)f_{0}(i-1)\Bigr{)}$$

$$f_3^{\text{ex}}(i^{(1)})$$

$$\displaystyle+g_{i}\Bigl{(}f_{0}(i)f_{2}(i-1)+f_{1}(i)f_{0}(i-1)+f_{2}(i)f_{1}(i-1)\Bigr{)}$$

$$H_{\mathrm{ch}}=\lambda_{\text{chirality}}\sum_{a}\sum_{i=2}^{N-1}\biggl{(}1-f_{a}(i^{(1)})\biggr{)}f_{a}^{\text{ex}}(i^{(1)})\,,$$

$$H_{\text{in}}=H_{\text{in}}^{(1)}+H_{\text{in}}^{(2)}+H_{\text{in}}^{(3)}+\dots$$

$$(j-i)=\{\begin{array}{ll}0\mathrm{mod}2& \mathrm{if}{\sum }_a\Delta n_a(i,j)=0\\ 1\mathrm{mod}2& \mathrm{if}{\sum }_a\Delta n_a(i,j)=(-1{)}^i\end{array}.$$

$$H^{(1)}=\sum _{i=1}^{N-4}\sum _{\begin{array}{c}j\ge i+5\\ j-i=1\mathrm{mod}2\end{array}}^N{h}_{i,j}^{(1)}$$

$$h^{(1)}_{i,j}=\mathrm{q}^{(1)}_{ij}\biggl{(}\epsilon^{(1)}_{ij}+\lambda_{1}(d(i,j)-1)+\sum_{r\in\mathcal{N}(j)}\lambda_{2}(2-d(i,r))+\sum_{m\in\mathcal{N}(i)}\lambda_{2}(2-d(m,j))\biggr{)}$$

$$H^{(2)}=H^{(2)}[1]+H^{(2)}[3,3]+H^{(2)}[3,5]+H^{(2)}[5,3].$$

$$H^{(2)}[1]=\sum _{i=1}^{N-4}\sum _{\begin{array}{c}j\ge i+4\\ j-i=0\mathrm{mod}2\end{array}}^N{h}_{ij}^{(2)}[1]$$

$$h_{ij}^{(2)}=\sum _{r\in \mathcal{N}(j)}{\text{q}}_{i,r}^{(1)}{ϵ}_{i,j}^{(2)}$$

$$H^{(2)}[3,5]=\sum _{i=1}^{N-4}\sum _{\begin{array}{c}j\ge i+4\\ j-i=0\mathrm{mod}2\end{array}}^N{h}_{ij}^{(2)}[3,5]$$

$$h_{ij}^{(2)}[3,5]=(1-\mathrm{q}^{(1)}_{i,r})(1-\mathrm{q}^{(1)}_{i,m})\mathrm{q}_{ij}^{(2)}[3,5]\biggl{(}\epsilon_{ij}^{(2)}+\lambda_{2}(d(i,j)-2)+\lambda_{3}(3-d(i,r))+\lambda_{5}(5-d(i,m))\biggr{)}$$

Full text

**Predicting the three-dimensional (3D) structure of a protein from its primary sequence of amino acids is known as the protein folding (PF) problem. Due to the central role of proteins’ 3D structures in chemistry, biology and medicine applications (e.g., in drug discovery) this subject has been intensively studied for over half a century Levitt and Warshel (1975); Zwanzig et al. (1992); Hinds and Levitt (1992); Shakhnovich (1995); Onuchic and Wolynes (2004); Lindorff-Larsen et al. (2011); Hyeon and Thirumalai (2011). Although classical algorithms provide practical solutions, sampling the conformation space of small proteins, they cannot tackle the intrinsic NP-hard complexity of the problem Unger and Moult (1993), even reduced to its simplest Hydrophobic-Polar model Berger and Leighton (1998). While fault-tolerant quantum computers are still beyond reach for state-of-the-art quantum technologies, there is evidence that quantum algorithms can be successfully used on Noisy Intermediate-Scale Quantum (NISQ) computers Moll et al. (2018); Preskill (2018) to accelerate energy optimization in frustrated systems Mandrà and Katzgraber (2018); King et al. (2017a, b); Streif and Leib (2019). In this work, we present a model Hamiltonian with $\mathcal{O}(N^{4})$ scaling and a corresponding quantum variational algorithm for the folding of a polymer chain with $N$ monomers on a tetrahedral lattice. The model reflects many physico-chemical properties of the protein, reducing the gap between coarse-grained representations and mere lattice models. We use a robust and versatile optimisation scheme, bringing together variational quantum algorithms specifically adapted to classical cost functions and evolutionary strategies (genetic algorithms), to simulate the folding of the 10 amino acid Angiotensin peptide on 22 qubits. The same method is also successfully applied to the study of the folding of a 7 amino acid neuropeptide using 9 qubits on an IBM Q 20-qubit quantum computer. Bringing together recent advances in building gate-based quantum computers with noise-tolerant hybrid quantum-classical algorithms, this work paves the way towards accessible and relevant scientific experiments on real quantum processors. **

The solution of Levinthal’s paradox Levinthal (1968) through a bias search in configuration space Zwanzig et al. (1992) demands a very fine description of the interactions in the cellular environment to correctly drive the search in the rugged energy landscape of a protein Onuchic et al. (1995); Onuchic and Wolynes (2004); Piana et al. (2014). GPU-assisted sampling methods of well parametrised coarse-grained models can provide useful insights regarding the protein’s native conformation, but folding a small protein in state-of-the-art simulations comes at a very high computational cost Duan (1998); Lindorff-Larsen et al. (2011). Protein lattice models reduce the conformational space and obviate the high computational cost of an off-lattice exhaustive sampling Levitt and Warshel (1975); Hinds and Levitt (1992). Quantum algorithms cannot ignore those simplifications because of curently available quantum resources. Perdomo-Ortiz et al. paved the way towards the construction of spin Hamiltonian to find the on-lattice heteropolymer’s low-energy conformations using quantum devices, but with unattainable high costs for NISQ computers Babbush et al. (2012); Perdomo Ortiz et al. (2008). Recently, the amount of resources needed for the simulation of polymer lattice models was reduced, however still maintaining an exponential cost in terms of number of qubits and gates needed Perdomo-Ortiz et al. (2012); Babej et al. (2018). These methods were used to fold a coarse-grained protein model with $6$ and $8$ amino acid sequences on a 2D and 3D lattice, respectively, using a quantum annealer Cai et al. (2014). These experiments required $81$ and $200$ qubits and led to a final population of $0.13\%$ and $0.024\%$ for the corresponding ground state structures, using divide and conquer strategies. More recently, Fingerhuth and coworkers proposed another approach based on the Quantum Approximate Optimization Algorithm (QAOA) Farhi and Harrow (2016) using a problem-specific alternating operator ansatz to model protein folding Fingerhuth et al. (2018). Employing the same model proposed in Babej et al. (2018) they succeeded in folding a $4$ amino acid protein model on a 2D square lattice.

In this work, we present a coarse-grained model for protein folding which is suited to the representation of branched heteropolymers comprised of $N$ monomers on a tetrahedral (or “diamond”) lattice. This choice is motivated by the chemical plausibility of the angles enforced by the lattice ($109.47^{o}$ for bond angles, $180^{o}$ or $60^{o}$ for dihedrals), which allows an all-atom description for a wide range of chemical and biological compounds. Given the modest resources of actual quantum devices, a two-centered coarse-grained description of amino acids (backbone and side chain) was used to mimic the protein sequence. Every monomer is depicted by one or multiple beads that can have a defined number of ‘color shades’ corresponding to different physical properties like hydrophobicity and charge.

The configuration qubits.

As for the previous models in literature, a polymer configuration is grown on the lattice by adding the different beads one after the other and encoding, in the qubit register the different “turn” $t_{i}$ that defines the position of the bead $i+1$ relatively to the previous bead $i$. Using a tetrahedral lattice, we distinguish two sets of nonequivalent lattice points $\mathcal{A}$ and $\mathcal{B}$ (see Fig. 1). At the $\mathcal{A}$ sites, the polymer can only grow along the directions $t_{i}\in\{0,1,2,3\}$ while at site $\mathcal{B}$ the possible directions are $t_{i}\in\{\bar{0},\bar{1},\bar{2},\bar{3}\}$. Along the sequence, the $\mathcal{A}$ and $\mathcal{B}$ sites are alternated so that we can use the convention that $\mathcal{A}$ (respectively $\mathcal{B}$) sites correspond to even (odd) $i$. Without loss of generality, the first two turns can be set to $t_{1}=\bar{1}$ and $t_{2}=0$ due to symmetry degeneracy. To encode the turns, we assign one qubit per axis $t_{i}=\mathrm{q}_{4i-3}\mathrm{q}_{4i-2}\mathrm{q}_{4i-1}\mathrm{q}_{4i}$ (Fig. 1(c)). Therefore, the total number of qubits required to encode a conformation $\bm{\mathrm{q}}_{\text{cf}}$ corresponds to $N_{\text{cf}}=4(N-3)$. If the monomers are described by more than one bead, the same formula holds by replacing $N$ with the total number of beads in the polymer. A denser encoding of the polymer chain using only $2(N-3)$ configuration qubits is presented in the SI.

The interaction qubits.

To describe the interactions, we introduce a new qubit register $\bm{\mathrm{q}}_{\text{in}}$, composed of $\mathrm{q}^{(l)}_{i,j}$ for each $l^{\mathrm{th}}$ nearest neighbour ($l$-NN) interaction on the lattice (see red and green dashed lines for $l=1$ and $l=2$ in Fig. 1(b)) between beads $i$ and $j$. The use of these registers will be explained in connection to the definition of the interaction energy terms. The number of qubits constituting the interaction register, $N_{\text{in}}$, is entirely determined by the skeleton of the polymer (i.e. including the side chains), regardless of the beads’ color, and scales as $\mathcal{O}(N^{2})$. Note that two $1$-NN beads occupy positions on different sub-lattices ($\mathcal{A}$ or $\mathcal{B}$). On the other hand, for $l>1$ all beads of both sub-lattices can potentially interact. Given a primary sequence, the pairwise interaction energies $\epsilon^{(l)}_{i,j}$ between the beads at distance $l$ can be arbitrarily defined to reproduce a fold of interest or it can be adapted from pre-existing models, like the one proposed by Miyazawa and Jernigan (MJ) for $1$-NN interactions Miyazawa and Jernigan (1996).

The Hamiltonian.

The next step defines the qubit Hamiltonian that describes the energy of a given fold defined by the sequence of beads (fixed) and the encoded turns. Penalty terms are applied when physical constraints are violated (e.g., when beads occupy the same position on the lattice), and physical interactions (attractive or repulsive in nature) are applied when two beads occupy neighbouring sites or are at distance $l>1$, where $l$ is the number of NN in real space. The different contributions to the polymer Hamiltonian are therefore (with $\bm{\mathrm{q}}=\{\bm{\mathrm{q}}_{\text{cf}},\bm{\mathrm{q}}_{\text{in}}$}),

[TABLE]

The definitions of the geometrical constraint ($H_{\text{gc}}$, which governs the growth of the primary sequence with no bifurcation) and the chirality constraint ($H_{\text{ch}}$, which enforces the correct stereochemistry of the side-chains if present) are given in SI.

The interaction energy terms.

For each bead $i$ along the sequence the distance to the other beads $j\neq i$ can uniquely be determined by the state of the $N_{\text{cf}}$ configuration qubits. To this end, for each pair of beads $(i,j)$ we introduce a four-dimensional vector (Eq. SI-13), the norm of which uniquely encodes they reciprocal distance $d(i,j)$. As an example, we consider the energy contributions for $1$-NN interactions. For each pair of beads $(i,j)$ an energy contribution of $\epsilon^{(l)}_{ij}$ is added to $H_{\rm{in}}^{i,j}$ when the distance $d(i,j)=l$. However, a contribution of the form $\epsilon^{(l)}_{ij}\,\delta(d(i,j)-l)$ cannot be efficiently implemented as a qubit string Hamiltonian (here $\delta(.)$ stands for the Dirac delta function). Using the set of contact qubits $\mathrm{q}^{(l)}_{i,j}$ we therefore define an energy term of the form $\mathrm{q}^{(l)}_{i,j}(\epsilon^{(l)}_{ij}+\lambda(d(i,j)-l))$ for each value of $l$ and $\lambda\gg\epsilon^{(l)}_{ij}$. This definition implies that the contribution $\epsilon^{(l)}_{ij}$ for the formation of the “interaction” $(i,j)$ at distance $l$ is only assigned when the contact qubit $\mathrm{q}^{(l)}_{i,j}=1$ and $d(i,j)=l$, simultaneously. For $\mathrm{q}^{(l)}_{i,j}=1$ and $d(i,j)\neq l$ the factor $\lambda$ adds a large positive energy contribution that overcomes the stabilizing energy $\epsilon^{(l)}_{ij}$ (the case of $d(i,j)<l$ is detailed in the SI).

Finally, in our model we prevent the simultaneous occupation of a single lattice site by two beads, as discussed in the SI. In a nutshell, we only prevent overlaps that occur in the vicinity of an interaction pair. If $\mathrm{q}^{(l)}_{i,j}=1$, we apply penalty functions so that $i$ and $j+1$ cannot overlap when $l=1$, for instance.

The folding algorithm.

The solution to the folding problem is the ground state of the Hamiltonian $H(\bm{\mathrm{q}})$ and therefore lies in the $2^{N_{\text{cf}}}$ dimensional space of the configuration qubits. To find this solution we prepare a variational circuit, comprising both the configurational and the interaction registers, which is composed by an initialization block with Hadamard gates and parametrized single qubit $R_{\text{Y}}$ gates followed by an entangling block and another set of single qubit rotations. We denote by $\bm{\theta}=(\bm{\theta}^{\text{cf}},\bm{\theta}^{\text{in}})$ the set of angles of size $2n$ where $n=N_{\text{cf}}+N_{\text{ct}}$ is the total number of qubits. Differently to the quantum mechanical case, for the solution of the ‘classical problem’ (e.g., folding) we do not need an estimate of the Hamiltonian expectation value, but we only require the sampling of the low energy tail of the energy distribution. Therefore, the optimization of the angles $\bm{\theta}$ is performed using a modified version of the Variational Quantum Eigensolver (VQE) Peruzzo et al. (2014); McClean et al. (2016) algorithm named Conditional Value-at-Risk (CVaR) VQE or simply CVaR-VQE Barkoutsos et al. (2019). Briefly, CVaR defines an objective function based on the average over the tail of a distribution delimited by a value $\alpha$ (see histogram in Fig. 2(a)) which is denoted CVaR${}_{\alpha}(\bm{\theta})=\langle\psi(\bm{\theta})|H(\bm{\mathrm{q}})|\psi(\bm{\theta})\rangle_{\alpha}$. Compared to conventional VQE, CVaR-VQE provides a drastic speed-up to the optimization of diagonal Hamiltonians as shown in Barkoutsos et al. (2019). The classical optimization of the gate parameters is performed using a Differential Evolution (DE) optimizer Storn and Price (1997), which mimics natural selection in the space of the angles $\bm{\theta}$. The optimisation procedure is summarized in Fig. 2(a). Note that at each step of the optimization, the wavefunctions $|\psi(\bm{\theta}^{p})\rangle$ corresponding to the different individuals $\bm{\theta}^{p}$ (Fig. 2(a)) are collapsed during measurement leading to binary strings, which are uniquely mapped to the corresponding configurations and energies. We denote by $\mathbb{P}_{\text{f}}(p)$ the probability for the individual $p$ to find the $\text{f}^{\text{th}}$ lowest energy fold at convergence.

Scaling.

We define the scaling of the algorithm as the number of terms (or Pauli strings), in the $n$-qubit Hamiltonian $H(\bm{\mathrm{q}})$ (see also Table I of SI).

[TABLE]

where $h_{\bm{\gamma}}$ real coefficients, $\mathrm{q}_{i}=(1-\sigma_{i}^{z})/2$ where $\sigma_{i}^{z}$ is $Z$ Pauli matrix with $\gamma_{i}\in\{0,1\}$, and $N_{\mathrm{t}}$ is the total number of terms. A thorough investigation of the scaling (see SI) reveals that the geometrical constraints imposed by the tetrahedral lattice give rise to all possible $2$-local terms within the $N_{\text{cf}}$ conformation qubits. Due to the coupling (entanglement) with the interaction qubits the Hamiltonian locality (i.e. the maximum number of Pauli operators different from the identity in $H(\bm{\mathrm{q}})$) is strictly $3$ for the $1$-NN interaction. Moreover, the scaling is bound by $N_{\mathrm{t}}\sim N_{\text{in}}{N_{\text{cf}}\choose{2}}=\mathcal{O}(N^{4})$ even for $l$-NN interactions, with $l\geq 1$. Fig. 1(d) and (e) respectively report the scaling of the proposed model and its qubits requirements.

Applications.

We first apply our quantum algorithm to the simulation of the folding of the 10 amino acid peptide Angiotensin. Using our coarse grained model on the tetrahedral lattice the simulation of this system would require 35 qubits, which is computationally intensive. We therefore introduced a denser encoding of the polymer configuration that requires only 2 qubits per turn $t_{i}=\mathrm{q}_{2i-1}\mathrm{q}_{2i}$, reducing the total number of qubits to 22. This variant generated $5$-local (instead of $3$-local) terms in the qubit Hamiltonian while keeping the total number of Pauli strings within an affordable range for small instances (see Fig. 1(d)). To further reduce the number of qubits we also integrate the side chains with the corresponding bead along the primary sequence and neglect interactions with $l>1$. Each bin of the histogram in Fig. 2(b) counts, over the population, the occurrence of $\mathbb{P}_{\text{f}}(p)$ for the minimum energy fold ($\mathrm{f}=0$) and the next 18 folds (histogram bars) given $n_{\mathrm{s}}=128$ (blue bars) and $n_{\mathrm{s}}=1024$ (orange bars) measurements of the wavefunctions during the minimization. More than $80\%$ of the individuals in the final population can generate the minimal conformations after $80$ generations (orange bars), which occurs with a probability $\max_{p}\mathbb{P}_{\text{f}}(p)=42.2\%$. The evolution of the percentage throughout the minimization can be found in the SI. By reducing the number of measurements to 128 shots, we obtained a broader spectrum of low energy conformations, which still includes the global minimum but with a lower probability. Among the low-energy conformations (with energies below 0), we can clearly identify the formation of an $\alpha$-helix and a $\beta$-sheet (conformations marked with a grey arrow in Fig. 2(b)). By tuning the interaction matrix (see Fig.S2 in SI), we can foster the formation of secondary structural elements. The 22-qubit Angiotensin system is still too large for encoding in state-of-the-art quantum hardware. To this end we used a smaller 7 amino acid neuropeptide with sequence APRLRFY (one-letter coding ) that can be mapped to 9 qubits. The corresponding CVaR-VQE circuit is shown in Fig. 3(a). As the entangling block we used a closed-loop of CNOT gates that fits the hardware connectivity of the 20-qubit IBM Q device Poughkeepsie (Fig. 3(b)). The mean CVaRα energy value of the population as a function of the number of generations shows a robust and smooth convergence towards the optimal fold (Fig. 3(d)). More importantly, the average probability (Fig. 3(e)) of the ground state averaged over the entire population, $\langle\mathbb{P}_{0}(p)\rangle$, increases monotonically reaching a final value larger than 20% and with $\max_{p}\mathbb{P}_{0}(p)$ peaking up at $33\%$ (see Fig. 3(c) and (e)).

Discussion and conclusions.

In this work, we introduced a quantum algorithm for the solution of the PF problem on a regular tetrahedral lattice. The model Hamiltonain describes a primary coarse-grained protein sequence where each beads represents an amino acid. Side chains can also be modeled by means of an additional bead linked to the main chain. The interaction between the amino acids (backbone and side chains) can be extended to $l$-NN (with $l>1$) along the lattice edges. This enables the modeling of sophisticated coarse-grained models accounting for Lennard-Jones and Coulombic like interactions. We show how the model can correctly reproduce secondary structure elements through the simple adaptation of the imputed contact map. The number of qubits scales quadratically with the number of amino acids , while the number of elements in the Hamltonian scales in $\mathcal{O}(N^{4})$. This implies the use of an unconventional treatment of the overlaps, which are avoided through the addition of penalty terms. Even though the PF problem is a classical optimization problem, the variational quantum algorithm used, CVaR-VQE, drastically reduces the number of measurements required to minimize the classical cost function (instead of the quantum mechanical average) and may lead to quantum advantage through the use of entanglement. The construction of specific mixing ansatz can drastically speed-up the search in the configuration space even when the ground state is not entangled but classical Hadfield et al. (2019); Fingerhuth et al. (2018). The direct connection between qubits and physical properties (configuration and contacts or interactions) allows a rationalisation of initialisation of the qubit states and their entanglement, beyond the simple scheme adopted in this preliminary investigation.

The locality of the Hamiltonian combined with the favourable scaling of the qubit resources and the circuit depth with the number of monomers, make our model the candidate of choice for the solution of the PF problem on NISQ devices and other quantum technologies. In this work we performed the simulation of the folding of the 10 amino acids protein Angiotensin using a realistic model for the noise of the one and two qubit gate operations. Furthermore we used a 20 qubit IBM Q processor to compute the folding of a 7 amino acid peptide on 9 qubits, which is to our knowledge the largest folding calculation on a NISQ device using a variational algorithm. The success of this calculation demonstrates the potential of our folding algorithm and opens up new interesting avenues for the use of quantum computers in the optimization of classical cost functions using the CVaR-VQE approach combined with a genetic algorithm for the selection of the best fitting variational parameters.

Author’s Contributions

A.R., S.W. and I.T. designed the project. A.R. and P.B. performed the experiments and the simulations. All authors contributed to the analysis of the results and to the writing of the manuscript.

Acknowledgements

I.T. and P.K.B. acknowledge financial support from the Swiss National Science Foundation (SNF) through the grant No. 200021-179312. A. R. is much obliged to Vladimir Nikolaevitch Smirnov who financially supported his work. All authors would like to acknowledge support from the IBM Q network and thank the qiskit development team for discussions regarding the development of the software.

IBM, IBM Q, Qiskit are trademarks of International Business Machines Corporation, registered in many jurisdictions worldwide. Other product or service names may be trademarks or service marks of IBM or other companies.

Materials and Methods

The noisy simulations were conducted with $\alpha=1\%$ (resp. $\alpha=0.1\%$) for the $128$ shots (resp. $1024$ shots) simulation with an “all-to-all“ entangling scheme on Qiskit Aleksandrowicz et al. (2019). All circuits were constructed with a VQE depth of $m=2$. Given a run with $n$ qubits, the size of the population for the evolutionary algorithm was set to $P=5mn$, a typical size according to the literature. The selection strategy of the DE algorithm is practically identical to the original “current-to-best/1/bin” Das et al. (2016).

Supporting Information for “Resource-Efficient Quantum Algorithm for Protein Folding”

I Lattice model

We call $i$ the backbone bead of index $i\in\{1,\dots,N\}$ in the primary sequence of the polymer and denote by $i^{(1)},i^{(2)},...,i^{(s)}$ the beads constituting its side chain. Without additional information, discussions about the locality of the Hamiltonian are made with respect to the sparser encoding.

I.1 Conformation encoding

Sparser Encoding

The polymer sequence is generated by specifying the series of turns on the lattice, starting from bead $i=1$ (see Fig. 1 of the main text). Each turn $t_{i}$ is encoded on $4$ qubits $\mathrm{q}_{4i-3}\mathrm{q}_{4i-2}\mathrm{q}_{4i-1}\mathrm{q}_{4i}$, each one representing a direction on the tetrahedral lattice (see Fig. 1 of the main text). One and only one qubit of the four will have value one (while the others will be set to zero). Accordingly, we encode the turns of the side chain in $i$ as $\mathrm{q}_{i^{(1)}}^{(0)}\mathrm{q}_{i^{(1)}}^{(1)}\mathrm{q}_{i^{(1)}}^{(2)}\mathrm{q}_{i^{(1)}}^{(3)},\mathrm{q}_{i^{(2)}}^{(0)}\mathrm{q}_{i^{(2)}}^{(1)}\mathrm{q}_{i^{(2)}}^{(2)}\mathrm{q}_{i^{(2)}}^{(3)}...$ where the upper index in parenthesis labels the 4 qubits describing the turn of side chain bead $i^{(k)}$ and the lower indices of the form $k^{(l)}$ label the bead $l$ along the side chain at the main chain position $k$. Without loss of generality, we choose the first two turns to be $0010$ and $0001$. A string of bits defining the entire conformation will have the general form of Eq. SI-1 where each ’residue’ (main and side chain beads) is embraced within square brackets and the side-chain qubits are within parentheses. Note that the first and the last bead of the main chain (terminal groups) have no side chains. This encoding requires $4(N-1)-8=4(N-3)$ qubits to totally define a conformation.

[TABLE]

Denser Encoding

The turn $t_{i}$ is encoded on two qubits $\mathrm{q}_{2i-1}\mathrm{q}_{2i}$. Accordingly, we will encode the turns in a given side chain by $\mathrm{q}_{i^{(1)}}^{(1)}\mathrm{q}_{i^{(1)}}^{(2)},\mathrm{q}_{i^{(2)}}^{(1)}\mathrm{q}_{i^{(2)}}^{(2)},...$ Without loss of generality we choose the first two turns to be $01$ and $00$. If the bead $2$ on the main chain does not bear a side chain, another qubit can be saved without breaking any symmetry ($\mathrm{q}_{6}=1$). The side chains can be encoded with the same convention as for the sparser encoding. A string of bits defining the entire conformation will have the general form of Eq. SI-2 where the side chain qubits are in parentheses. This encoding requires $2(N-1)-4=2(N-3)$ qubits to totally define a conformation.

[TABLE]

Turn indicator

It is convenient to introduce an indicator for the axis chosen at turn $t_{i}$. We will denote $f_{a}(\mathrm{q}_{2i-1},\mathrm{q}_{2i})=f_{a}(i)$ the function that returns $1$ if the axis $a\in\{0,1,2,3\}$ (see main text, Fig. 1) is chosen at turn $t_{i}$. For the denser encoding this function is given by

[TABLE]

For the sparser encoding, the expression of $f_{a}(i)$ is trivial:

[TABLE]

Distances

The shortest spatial distance between two beads is measured as a function of the number of turns separating them along the main chain (see below for the case where the distance is computed between two beads of the side chains). We define as $n_{a}(i,j)$ (resp. $n_{\bar{a}}(i,j)$) the number of occurrence of $a\in\{0,1,2,3\}$ (reps. $\bar{a}\in\{\bar{0},\bar{1},\bar{2},\bar{3}\}$) along the sequence separating $i$ and $j$ with $j>i$. All intermediate beads will be labeled by $k=i,\dots,j-1$. We assume that the polymer starts with a bead on the sub-lattice $\mathcal{B}$. Due to the alternation of the two sub-lattices, all sites of the sub-lattice $\mathcal{A}$ ($\mathcal{B}$) will be labelled with even (odd) integers. Let’s denote $\Delta n_{a}(i,j)=n_{a}(i,j)-n_{\bar{a}}(i,j)$ with $n_{a}(i,j)$ (resp. $n_{\bar{a}}(i,j)$) the number of occurrence of $a$ (reps. $\bar{a}$) in the sequence between bead $i$ and $j$. We then have

[TABLE]

The distance between beads that are placed on side chains, say $i^{(s)}$ ($s^{th}$ bead of side chain at position $i$ of the main chain) and $j^{(p)}$ ($p^{th}$ bead of side chain at position $j$) with $j>i$ is obtained by first considering the distance between $i$ and $j$ and then add or remove the contributions of side chain qubits. Because $j>i$, the side chain qubits of $j$ (resp. $i$) are taken into account with a plus (resp. minus) sign. Moreover, the direction is defined by the parity of $i$ and $j$. The generalized expression for $\Delta n_{a}$ is therefore

[TABLE]

For the calculation of the actual distance, is is convenient to introduce the four dimensional vector $\mathbf{x}(i,j)$ defined as

[TABLE]

such that $d(i,j)=||\mathbf{x}(i,j)||^{2}_{2}=\sum_{a}\Delta n_{a}(i,j)^{2}$. It is crucial to note that on the tetrahedral lattice there is a bijective map between the values of $d(i,j)$ and the Euclidean (through-space) distances $r_{ij}$ in lattice bond units,

[TABLE]

I.2 Construction of local Hamiltonians

Sparser encoding Hamiltonian

For the sparser encoding, we need to impose that one and only one of the four qubits that define a turn is equal to one. On a quantum circuit, this can be achieved easily by using a valid initialization of the qubits together with gate operations that conserved the number of 1’s in each set of 4 qubits that encodes a turn. When this is not possible, a penalty function as show in Eq SI-15 with a large positive $\lambda$ can be used to impose this constraint.

[TABLE]

Growth constraint.

Several constraints are needed to prevent the growth of the chain towards unphysical geometries (e.g., to prevent that the chain (main or side) at side $i$ folds back into itself). To this end, we compute function $T(i,j)$ defined as

[TABLE]

for each pair of beads $i$ and $j$. $T(i,j)$ returns a 1 if and only if the turns $t_{i}$ and $t_{j}$ are along the same axis ($a$ and $\bar{a}$, respectively). Note that $T(i,j)$ is composed of $2$-local terms for the sparse encoding. Firstly, we need to eliminate sequences where the same axis ($a$ and $\bar{a}$) is chosen twice in a row (e.g. $2\bar{2}$) since this will give rise to a chain folding back into itself. To this end we apply the following penalty term

[TABLE]

with large positive $\lambda_{\text{back}}$. Note that we can easily control the number of 2-local terms appearing in the sum (linear number of two local terms). In the general case, for a degree of branching $s>1$, similar terms need to be added to prevent the overlap of two consecutive bonds within the side chains.

Chirality constraints

Natural polymers have a well defined chirality that has to be imposed in our model. In proteins, the position of the side chains at the insertion point with the main chain determines the chirality of each residue. The position of the first side chain bead $i^{(1)}$ on the main bead $i$ is imposed by the choice of the (main chain) turns $t_{i-1}$ and $t_{i}$.

To enforce the correct chirality, we add a constraint $H_{\mathrm{ch}}$ to the Hamiltonian. The required (expected) chirality at $i^{(1)}$ is encoded in the function $f_{a}^{\text{ex}}(i^{(1)})$ which is a function of $f_{a}(i-1)$ and $f_{a}(i)$ defined in Eq.SI-6.

Analyzing all possible cases, we can construct the truth-table given in Table S1, which uniquely define $f_{a}^{\text{ex}}(i^{(1)})$ for $a\in\{0,1,2,3\}$. To this end, we first need to define an indicator for the parity, $g_{i}$, which is 0 if $i$ is even and 1 otherwise. We then write $g_{i}=\frac{1-(-1)^{i}}{2}$. Using the information in Table S1 and the definition of $g_{i}$, we obtain

[TABLE]

The penalty term to impose the right chirality, $H_{\mathrm{ch}}$, will then the penalize the structures for which $f_{a}(i^{(1)})$ differs from $f_{a}^{\text{ex}}(i^{(1)})$,

[TABLE]

for large and positive values of $\lambda_{\text{chirality}}$. Note that this term actually only contains linear number of 3-local terms for the denser encoding (5-local for the sparser encoding).

I.3 Construction of $H_{\text{in}}$

Definitions.

We can decompose the interaction Hamiltonian into different Hamiltonian terms $H^{(l)}$ corresponding to $l^{t}h$ nearest neighbour interactions,

[TABLE]

To keep the number of terms in $H_{\text{in}}$ as low as possible we will need to exploit properties of the tetrahedral lattice. In this section, we only consider interactions between beads $i$ and $j$ with $j>i+2$. Indeed, interactions between beads that are nearest or second nearest neighbour in the primary sequence are present in every conformations so considering them will not help discriminating the folds. Note also that beads that are separated by less than 5 bonds cannot be nearest neighbour (1-NN) on the lattice. These beads can indifferently represent main chain beads as well as side chain beads since only the absolute distance measure in bond units along the polymer chain matters (e.g. $j^{(1)}-i\equiv j-i+1$ or $j^{(1)}-i^{(1)}\equiv j-i+2$). First, we observe that the parity of $j-i$ corresponds to the parity of $d(i,j)$. Then, since we alternatively choose among positive and negative directions to encode a conformation, we have

[TABLE]

In addition, we introduce a notation for the set of first nearest neighbor beads to a given bead $i$, which will be denoted $\mathcal{N}(i)$. For example, $\mathcal{N}(i)=\{i-1,i+1,i^{(1)}\}$ if $i$ is not at the end of the polymer chain or $\mathcal{N}(i^{(1)})=\{i\}$ when limited to one side chain bead per monomer. Finally, we define a new set of two-index quantities $\mathrm{q}^{(l)}_{ij},\mathrm{q}^{(l)}_{i^{(1)}j},\mathrm{q}^{(l)}_{i,j^{(1)},}\mathrm{q}^{(l)}_{i^{(1)}j^{(1)}}$ that encode the presence of an interaction of order $l$ between beads $i$ and $j$ (i.e, $i$ and $j$ are l-NN).

First nearest neighbour interactions: $H^{(1)}$.

The quantity $\mathrm{q}^{(1)}_{ij}$ describes a contact between a pair of beads $i$ and $j$ that belong to different lattices and that can be any far apart along the polymer chain. If this contact occurs, the energy will be modified by the contact energy $\epsilon^{(1)}_{ij}$. The corresponding term in the Hamiltonian is then $\mathrm{q}^{(1)}_{ij}\epsilon^{(1)}_{ij}$. However, since we keep the configuration and the interaction qubit registers independent from each other, before assigning the energy contribution $\epsilon^{(1)}_{ij}$ we need to make sure that the two beads are indeed at a distance of $1$. The correct energy term is therefore $\mathrm{q}^{(l)}_{ij}=1$ and reads $\mathrm{q}^{(1)}_{ij}\lambda_{1}(d(i,j)-1)$, which implies that we apply either the stabilization energy $\epsilon^{(1)}_{ij}$ when both $\mathrm{q}^{(1)}_{ij}=1$ and $d(i,j)=1$, or a penalty contribution weighted by the positive factor $\lambda_{1}$. Here $d(i,j)$ cannot be equal to 0 because the two beads must be on different lattices. Note that if $r\in\mathcal{N}(j)$ and $i$ and $j$ are first nearest neighbours, $r$ can either overlap with $i$ or be at a distance of $2$ from $i$. In order to avoid the overlap of different beads, we penalize the occurrence of the first option by applying a penalty function. This constraint can be written as $\mathrm{q}^{(1)}_{ij}\lambda_{2}(2-d(i,j))$ with $\lambda_{2}$ being a large positive value; it will prevent the formation of a local overlaps between polymer beads in the vicinity of a nearest neighbor contact. Summing over all the possible contacts between backbone beads, we obtain the general form of the first nearest neighbour interaction Hamiltonian

[TABLE]

with

[TABLE]

The choice of the ratio $\lambda_{1}/\lambda_{2}$ has to be made wisely. In fact, we need that the term within parentheses in Eq. (SI-30) is always positive when $i$ and $j$ are not in contact. Since $|\mathcal{N}(i)|\leq 3$ and $-(j-i+1)<d(i,r)-d(i,j)<(j-i+1)$ when $r\in\mathcal{N}(j)$, we need $\lambda_{1}>6(j-i+1)\lambda_{2}+\epsilon^{(1)}_{ij}$. The penalty values $\lambda_{1}$ and $\lambda_{2}$ therefore depend on $i$ and $j$.

Second nearest neighbour interactions: $H^{(2)}$.

In order to take into account second nearest neighbours interactions between beads $i$ and $j$ (which by definition occupy the same sub-lattice), we need to introduce a new set of interaction qubits: $\mathrm{q}^{(2)}_{ij}$. In this case, there are multiple possible configurations that lead to 2-NN interactions, so multiple interaction qubits are needed for each pair $(i,j)$. We deal with these different cases by splitting the Hamiltonian as follows

[TABLE]

$H^{(2)}[1]$ deals with the case in which a given configuration leads to a direct contact between beads $r$ and $i$ or $m$ and $i$ where $(r,m)\in\mathcal{N}(j)$ (see Fig S1),

[TABLE]

where

[TABLE]

The remaining configurations can be classified into the following classes depending of the distances between bead $i$ and the beads $r$ and $m$, with $r,m\in\mathcal{N}(j)$ (see Fig S1):

• •

$[3,3]:d(i,r)=3$ and $d(i,m)=3$

• •

$[3,5]:d(i,r)=3$ and $d(i,m)=5$

• •

$[5,3]:d(i,r)=5$ and $d(i,m)=3$.

We thus introduce three second contact qubits $\mathrm{q}_{ij}^{(2)}[3,3],\mathrm{q}_{ij}^{(2)}[3,5],\mathrm{q}_{ij}^{(2)}[5,3]$ and use the same method to build the Hamiltonian. For instance,

[TABLE]

with

[TABLE]

Note that we also add a check for $\mathrm{q}^{(1)}_{i,r}\neq 1$ and $\mathrm{q}^{(1)}_{i,m}\neq 1$ to avoid double counting of the interaction energy $\epsilon_{i,j}^{(2)}$. Also in this case, we choose $\lambda_{2},\lambda_{3},\lambda_{5}$ so that the term in the parenthesis is always positive if $d(i,j)\neq 2$. Additional qubits are necessary when dealing with side chains. The nature of these additional terms can be easily elaborated once the gist of the method is understood. To take into account $l$-NN interactions with $l>2$ we follow the same method of enumerating all the possible configurations giving rise to the specific interactions. The number of required qubits will evolve exponentially with $l$ but still quadratically with $N$ if truncated to a cutoff maximum value. The interaction energy can also depend on the configuration of interaction.

II Contact Map

Contact map for the stabilization of two different secondary structure elements: $\beta$-sheet (left) versus $\alpha$-helix (right). Simulations performed with these contact maps reproduce the correct minimal energy structures (see corresponding cartoons). The $\alpha$-helix is stabilized by a second line of contacts parallel to the main diagonal, while the (anti-parallel) $\beta$-sheet forms a series of contacts along the counter-diagonal.

III Convergence of CVaR VQE for noisy simulation

Simulation of the folding of the peptide Angiotensin using a realistic noise model for hardware.