**Ground-State Protein Folding Using Variational Quantum Eigensolver (VQE)**

This example shows an efficient method for using qubits to encode a protein fold on a 3-D tetrahedral lattice [1], [2]. The ground-state is found through a simulated variational quantum eigensolver (VQE) routine. The VQE algorithm uses classical optimization to improve the initial guess of the ground state, and then a quantum computer calculates the expectation value. The final circuit from the simulation is run on a real QPU for comparison.

The protein is a neuropeptide with seven amino acids, APRLRFY, pictured below. The example assumes a coarse-grained protein model, where "beads" representing amino acids can traverse the lattice and interact with each other. Each bond between amino acids can be in one of four directions, corresponding to corners of a tetrahedron. These four turns can be represented by two qubits.

**Model Protein Fold with Qubits**

The purpose of the embedding methods of [1] and [2] is to model the physical pairwise interaction energy of the protein beads while also imposing energy penalties on unrealistic configurations (such as overlaps or beads on the same lattice site). Therefore, the example is configurable in that you can modify the amino acid letters to test the effects on folding. However, the example code assumes a protein with seven amino acids, so you cannot change the length of the protein string.

**Configuration Qubits**

A protein of `N`

beads (here `N=7`

) can make `N-1`

turns, therefore `2*(N-1)=12`

bits are required for the APRLRFY protein bonds. However, without loss of generality, the first two turns can be fixed to `01`

and `00`

. One of the bits in the third turn is fixed due to other symmetry considerations. The `turn2qubit`

mapping depicts the 12 bits, denoting the value of the 5 that are fixed, and the 7 that are variable which will be represented by qubits. For details, see the dense encoding scheme in [2].

hyperParams.protein = 'APRLRFY'; hyperParams.turn2qubit = '0100q1qqqqqq'; hyperParams.numQubitsConfig = sum(hyperParams.turn2qubit=='q');

**Interaction Qubits**

The methods of [1] and [2] can consider interactions between an arbitrary number of nearest-neighbors (NN). This example only considers 1-NN interaction terms, which are only possible between some beads due to the structure of the lattice. Beads separated by less than 5 bonds cannot be 1-NN (Eq. 29 of [2]). There are two interaction pairs in the length-7 protein, between beads 1 and 6, and between beads 2 and 7. Two interaction bits are used to specify whether there is any interaction between the beads in each of these pairs.

hyperParams.numQubitsInteraction = 2;

A matrix of contact energy values is used as a look-up table for pairwise interaction between different amino acids. For this example, the values are chosen at random.

function mat = buildMJInteractions(protein) % Create the MJ interaction energy matrix for the protein, specified with % 1-letter codes N = length(protein); mat = zeros(N,N); rng(29507) MJ = rand(20)*-6; MJ = triu(MJ) + triu(MJ, 1)'; acids = ["C","M","F","I","L","V","W","Y","A","G","T","S","N","Q","D","E","H","R","K","P"]; acid2idx = dictionary(acids, 1:20); for i = 1:N for j=1:N mat(i,j) = MJ(acid2idx(protein(i)), acid2idx(protein(j))); end end end

Call the `buildMJInteractions`

function to find the interaction energies for the protein.

hyperParams.interactionEnergy = buildMJInteractions(hyperParams.protein);

**Write Function to Calculate Energy of Folds**

The `exactHamiltonian`

function calculates the energy for each possible fold of the protein. Only 1-NN interactions are considered, and energy penalties are imposed for unrealistic configurations.

function energies = exactHamiltonian(bitstrings,hyperParams) % Compute the Hamiltonian for each bit string (i.e., the energy for each fold). % See [2] for details. This does not consider the Hch constraint from % side-chains and the interaction term is only 1-nearest-neighbor (1-NN). % H = Hgc + Hin_1 lambdaDis = 720; % Penalty for interaction distance lambdaLoc = 20; % Penalty for interaction location lambdaBack = 50; % Penalty for unphysical geometry energies = zeros(size(bitstrings,1),1); numBeads = length(hyperParams.protein); for idx = 1:length(energies) bitstring = bitstrings(idx,:); config = hyperParams.turn2qubit; config(config=='q') = bitstring(1:hyperParams.numQubitsConfig); turns = bin2dec(reshape(config,2,numBeads-1)'); %% Geometric Hamiltonian Hgc % Count number of adjacent turns which are equal and impose a penalty for each energies(idx) = lambdaBack*sum(turns(1:end-1) == turns(2:end)); %% 1-NN Interaction Hamiltonian Hin currInteractionQubit = hyperParams.numQubitsConfig; for i=1:(numBeads-4) for j=(i+5):2:numBeads currInteractionQubit = currInteractionQubit+1; if bitstring(currInteractionQubit)=='0' continue; end % Add the interaction energy energies(idx) = energies(idx) + hyperParams.interactionEnergy(i,j); % Compute distances between interacting beads deltaN_ij = zeros(1,4); deltaN_ir = zeros(1,4); deltaN_mj = zeros(1,4); for k=0:3 deltaN_ij(k+1) = (-1).^(i:j-1)*(turns(i:j-1) == k); deltaN_ir(k+1) = (-1).^(i:j-2)*(turns(i:j-2) == k); deltaN_mj(k+1) = (-1).^(i+1:j-1)*(turns(i+1:j-1) == k); end d_ij = norm(deltaN_ij)^2; d_ir = norm(deltaN_ir)^2; d_mj = norm(deltaN_mj)^2; % Add penalty for distance not equal to 1 energies(idx) = energies(idx) + lambdaDis*(d_ij-1); % Add penalty for unphysical nearest neighbour collisons energies(idx) = energies(idx) + lambdaLoc*(2-d_ir); energies(idx) = energies(idx) + lambdaLoc*(2-d_mj); if i-1 >= 1 for k=0:3 deltaN_mj(k+1) = (-1).^(i-1:j-1)*(turns(i-1:j-1) == k); end d_mj = norm(deltaN_mj)^2; energies(idx) = energies(idx) + lambdaLoc*(2-d_mj); end if j+1 <= numBeads for k=0:3 deltaN_ir(k+1) = (-1).^(i:j)*(turns(i:j) == k); end d_ir = norm(deltaN_ir)^2; energies(idx) = energies(idx) + lambdaLoc*(2-d_ir); end end end end end

**Compute Minimum Energy for All Folds**

Nine qubits are required to fully model a fold (7 configuration qubits and 2 interaction qubits). Use the `exactHamiltonian`

function directly with all possible bit strings (each representing a possible fold), to find the minimum energy fold. You can later compare this value to the result of the quantum-based optimization. This exhaustive search is realistic for small proteins, but is not practical for larger proteins since the number of combinations of folds grows quickly.

hyperParams.numQubitsTotal = hyperParams.numQubitsConfig + hyperParams.numQubitsInteraction; allFolds = dec2bin(0:2^hyperParams.numQubitsTotal-1,hyperParams.numQubitsTotal); allEnergies = exactHamiltonian(allFolds,hyperParams); hyperParams.GroundState.Energy = min(allEnergies); hyperParams.GroundState.Index = find(allEnergies == hyperParams.GroundState.Energy); allFolds(hyperParams.GroundState.Index,:)

`ans = `*2×9 char array*
'101000110'
'101001010'

There are two lowest-energy folds because the last amino acid can occupy either position in the lattice. Find the energy of each of the identified folds.

allEnergies(hyperParams.GroundState.Index)

`ans = `*2×1*
-4.7973
-4.7973

This minimum-energy value matches the interaction energy between beads 1 and 6.

**Write CVaR-VQE Objective Function**

Next, use a quantum circuit to translate this optimization problem over a set of bit strings into a different optimization problem, this one over a set of angles between `-pi`

and `pi`

. These angles have no connection to the angles between the protein beads. The objective function accepts the set of angles for use inside of the quantum circuit, and the circuit returns various bit strings with different probabilities. Call the `exactHamiltonian`

function on the bit strings returned by the quantum circuit, and use a weighted mean of the smaller energies computed as the value of the objective function.

Specifically, after the energy is computed for each observed fold, the associated probabilities are sorted by energy. The objective function returns an expectation energy computed from the tail end of the probability distribution, cutoff by an alpha parameter. This expectation energy is a conditional value at risk (CVaR). An alpha value of 0.05 was used experimentally in [1], but for a noise-free simulation `ProteinVQEObjective`

uses a smaller cutoff value of 0.025.

function [energy, maxProbFold] = ProteinVQEObjective(parameters,hyperParams) % Construct and simulate the variational circuit ansatz = ProteinConfigAnsatz(parameters); qState = simulate(ansatz); allProbs = (qState.Amplitudes).^2; % There are 10 qubits in the circuit, but only these 9 are used to define a fold foldQubits = [1:8 10]; % Get the most probable fold % Only compute this if second output maxProbFold is requested if nargout > 1 [~,idx] = max(allProbs); maxProbKet = char(qState.BasisStates(idx)); maxProbFold = maxProbKet(foldQubits); end % Sample, and query/get the states and probabilities of the fold qubits qMeasurement = randsample(qState, hyperParams.numShots); [states, probs] = querystates(qMeasurement, foldQubits); % Sort the probabilities by the energy [energies,sort_idx] = sort(exactHamiltonian(char(states), hyperParams)); probs = probs(sort_idx); % Compute CVaR over the low energy tail of the energy distribution, % delimited by a cutoff parameter alpha. alpha = .025; cut_idx = nnz(cumsum(probs) < alpha); cvar_probs = probs(1:cut_idx); cvar_probs(end+1) = alpha - sum(cvar_probs); % Compute expectation energy as the sum of cutoff state energies weighted by their probability energy = dot(cvar_probs, energies(1:cut_idx+1))/alpha; end

Define the number of shots to use in `ProteinVQEObjective`

, and create a function handle that passes in all of the parameter values to `ProteinVQEObjective`

.

hyperParams.numShots = 1024; objFcn = @(theta) ProteinVQEObjective(theta,hyperParams);

### Create Circuit Ansatz

Create and plot the variational circuit ansatz for this protein with random angles. Qubits 1-7 represent the configuration, and the other three qubits are for the interaction (including a helper qubit that does not get measured). This circuit can be found as Figure 3 in [1].

function ansatz = ProteinConfigAnsatz(parameters) % Create the circuit ansatz for a 7 amino acid neuropeptide (10 qubit circuit). parameters = reshape(parameters, [2 9]); gates = [hGate([1:7 9 10]) ryGate([1:7 9 10], parameters(1,:)) cxGate(1:4, 2:5) cxGate(5,10) cxGate(10,9) cxGate(9,8) cxGate(8,9) cxGate(9,8) cxGate([8 7 6], [7 6 1]) ryGate([1:8 10], parameters(2,:)) ]; ansatz = quantumCircuit(gates); end

Call `ProteinConfigAnsatz`

with random angles to construct the circuit. Plot the circuit to view the qubits and gates.

ansatz = ProteinConfigAnsatz(rand(2,9)); plot(ansatz)

**Simulate Iterations of CVaR-VQE**

There are two layers of RY rotation gates over the configuration and interaction qubits. These rotation angles are learned parameters in the variational circuit. Use `surrogateopt`

(Global Optimization Toolbox) with the objective function to find the set of RY rotation angles that produce the minimum expectation energy. The original paper [1] uses a genetic algorithm to find the rotation angles. `surrogateopt`

converges to similar values, but because it uses intervals internally, the convergence behavior is comparatively different.

Define the number of angles and set options for the maximum number of function evaluations, plotting function, and initial points. Then call `surrogateopt`

with the objective function, upper and lower bounds for the angles, and options.

numAngles = 2*hyperParams.numQubitsTotal; rng default options = optimoptions("surrogateopt",... "MaxFunctionEvaluations",10, ... "PlotFcn","optimplotfval",... "InitialPoints",pi*ones(numAngles,1)); lb = repmat(-pi,numAngles,1); ub = repmat(pi,numAngles,1); [angles,minEnergy] = surrogateopt(objFcn,lb,ub,[],[],[],[],[],options);

surrogateopt stopped because it exceeded the function evaluation limit set by 'options.MaxFunctionEvaluations'.

Call the objective function using the optimized angles to find the state bit string with highest probability as well as its associated energy.

[groundStateEnergy,groundStateFold] = ProteinVQEObjective(angles,hyperParams)

groundStateEnergy = -4.7973

groundStateFold = '101000110'

This ground-state fold is one of the two folds obtained earlier with the minimum energy calculation.

allFolds(allEnergies==minEnergy,:)

`ans = `*2×9 char array*
'101000110'
'101001010'

Visualize the ground-state fold using the `plotProtein`

function, which enables you to plot the protein structure from the qubit values.

function plotProtein(bitstring,hyperParams) % Plot protein structure from the bitstring % The input bitstring is expected to be of length 9, with the first 7 bits % specifying turns in direction in the structure, and the last 2 bits % specifying interactions between beads. % Number of beads N = length(hyperParams.protein); % Construct 3D coordinates representing the 4 corners of a tetrahedron % centered in 0. These represent the 4 directions each new bead might % be added on in a tetrahedral grid. turn2bead = ones(4,3); turn2bead(2:4,:) = -1+2*eye(3); % Construct complete bitstring by inserting input bitstring into % the complete string mask. completeBitstring = hyperParams.turn2qubit; completeBitstring(completeBitstring=='q') = bitstring(1:hyperParams.numQubitsConfig); % Each pair of bits is converted into a number between 0 and 3 turns = bin2dec(reshape(completeBitstring,2,[])'); % Change direction in which we follow the tetrahedral grid in each bead % on the line. signs = (-1).^(0:N-1)'; % Compute placements of each bead. beads = cumsum(signs.*[zeros(1,3);turn2bead(turns+1,:)]); % Plot the beads connecting by lines, and add a text label to each figure plot3(beads(:,1),beads(:,2),beads(:,3),'.-','LineWidth',2,'MarkerSize',80,'SeriesIndex',1) axis off viewDir = -[6 4 1]; view(viewDir); % Plot text, ensuring it is in front of the lines beadsText = beads + 0.03*viewDir; text(beadsText(:,1),beadsText(:,2),beadsText(:,3),hyperParams.protein', ... 'FontWeight', 'bold', 'HorizontalAlignment','center') % Whether there are interactions between any pair of beads is determined % by the additional bits in the string. interactions = []; currInteractionQubit = hyperParams.numQubitsConfig+1; for i=1:(N-5) for j=(i+5):2:N if bitstring(currInteractionQubit) == '1' interactions = [interactions;beads(i,:);beads(j,:);nan*ones(1,3)]; %#ok<AGROW> end currInteractionQubit = currInteractionQubit+1; end end if ~isempty(interactions) hold on plot3(interactions(:,1),interactions(:,2),interactions(:,3),'k--','LineWidth',2) hold off legend([hyperParams.protein+" Protein Structure";"Interactions"], "Location","southoutside") else legend(hyperParams.protein+" Protein Structure", "Location","southoutside") end end

The protein has two possible interactions considered in the 1-NN model. The fold with the lowest energy only has one of these interactions. Use `plotProtein`

to visualize the lowest energy fold.

plotProtein(groundStateFold,hyperParams)

Next, construct the quantum circuit using the optimized angles and simulate the circuit to see the expected probability distribution over states. Use a threshold to filter out states with probabilities less than 2%. The ground-state fold `101000110`

appears again as the state with the highest probability.

optimized_circuit = ProteinConfigAnsatz(angles); sv = simulate(optimized_circuit); histogram(sv,[1:8 10],Threshold=0.02)

**Run Final Iteration on QPU**

Connect to the IonQ device using `quantum.backend.QuantumDeviceAWS`

. Specify the region of the device and path to a bucket to store results.

reg = "us-east-1"; bucketPath = "s3://amazon-braket-mathworks/doc-examples"; device = quantum.backend.QuantumDeviceAWS("Harmony",S3Path=bucketPath,Region=reg)

device = QuantumDeviceAWS with properties: Name: "Harmony" DeviceARN: "arn:aws:braket:us-east-1::device/qpu/ionq/Harmony" Region: "us-east-1" S3Path: "s3://amazon-braket-mathworks/doc-examples"

Create a task to run the circuit ansatz with optimized angles on the QPU. Specify the number of shots as 500.

task = run(optimized_circuit,device,NumShots=500); wait(task)

Fetch the results and plot a histogram of the states for all but the ninth qubit. Specify a lower threshold of 2% to filter out unlikely states.

results = fetchOutput(task); histogram(results,[1:8 10],Threshold=0.02)

The ground-state fold `101000110`

appears again as the state with the highest probability, so the QPU results agree with the local simulation of the circuit.

### References

[1] Robert, Anton, Panagiotis Kl. Barkoutsos, Stefan Woerner, and Ivano Tavernelli. “Resource-Efficient Quantum Algorithm for Protein Folding.” *Npj Quantum Information* 7, no. 1 (February 17, 2021): 38. https://doi.org/10.1038/s41534-021-00368-4.

[2] Robert, Anton, Panagiotis Kl. Barkoutsos, Stefan Woerner, and Ivano Tavernelli. “Supplementary Information for 'Resource-Efficient Quantum Algorithm for Protein Folding'" *Npj Quantum Information* 7, no. 1 (February 17, 2021): 38. https://doi.org/10.1038/s41534-021-00368-4.