Inverse Protein Folding Problem via Quadratic Programming

0.4 MB·

by  Andrii Riazanov

Inverse Protein Folding Problem
via Quadratic Programming
Andrii Riazanov Mikhail Karasikov
Skolkovo Institute of Science and Technology Moscow Institute of Physics and Technology
Moscow Institute of Physics and Technology Skolkovo Institute of Science and Technology
[email protected] [email protected]
Sergei Grudinin
7
1 University of Grenoble Alpes, LJK
0
CNRS, LJK
2
Inria
n
a [email protected]
J

3

]
M Abstract design (CPD), which targets at choosing a sequence of
amino acids to perform a particular task. As there are
Q
This paper presents a method of reconstruction a pri- 20 possible amino acids for each position in the protein
o. mary structure of a protein that folds into a given geo- chain, the variety of sequences for a specific 3D shape
i metricalshape. Thismethodpredictstheprimarystruc- grows exponentially. In CPD, the inverse folding prob-
b
ture of a protein and restores its linear sequence of lemisformulatedasanoptimizationproblem,aimedat
-
q amino acids in the polypeptide chain using the tertiary the minimization of the energy of the protein.
[ structure of a molecule. Unknown amino acids are de- The inverse folding problem is usually solved in-
1 termined according to the principle of energy minimiza- cluding the fact that each residue can adopt multiple
v tion. This study represents inverse folding problem as conformational states, which are most often grouped
3 a quadratic optimization problem and uses different re- into a finite number of conformational isomers (ro-
7
laxationtechniquestoreduceittotheproblemofconvex tamers). Two approximations are common for this
6
optimizations. Computational experiment compares the problem. First, it is assumed that the protein’s back-
0
0 quality of these approaches on real protein structures. bone is fixed. Second, the continuous domain of avail-
. able conformations for the side-chains is approximated
1
0 using a finite set of discrete conformations, defined by
7 1. Introduction inner dihedral angles. The CPD is then formulated
1 as the problem of determining the conformation cor-
:
v Aproteinisasequenceofbasicmoleculesofamino respondingtotheminimumenergyfromthatfiniteset.
i acids. Proteins are essential compounds of all living A variety of algorithms was developed to solve
X
organisms, they are involved in almost all structural, the inverse folding problem. In [1], the problem is
r
a catalytic, sensory, and regulatory functions. Proper- expressed as a Cost Function Network (CFN), also
ties and functions of proteins mostly depend on their known as the Weighted Constraint Satisfaction Prob-
structure. Protein engineering is aimed at studying of lem (WCSP). The interaction energies are represented
protein structures and it finds applications in a num- as cost functions, and the unknown variables are the
ber of areas. For example, the design of proteins with sets of pairs (amino acid, conformation). Then, the
desired shapes and properties is required in such areas problem is solved using the dead-end elimination/A∗
asmedicine,biotechnology,syntheticbiology,nanotech- algorithm (DEE/A∗). More precisely, first, the DEE
nologies, etc. algorithm reduces the computational cost of the prob-
While computational protein folding predicts the lem, excluding from the consideration conformations
structureofagivensequence,thegoalofcomputational that are unlikely to be optimal. Then, the A∗ algo-
inverse protein folding problemistofindaminoacidse- rithmallowstoexpandasequence-conformationtree,so
quences that would fold best into a given 3D shape or that sequence-conformations are extracted and sorted
scaffold. It can be referred as computational protein on the basis of their energy values. The CFN prob-
1
lem is solved using the toulbar2 solver. The investiga- [14, 15], which describe only interactions between all
tion [1] also compares the CFN approach with several possible pairs of amino acids, but not rotamers.
other optimization models. CPD is expressed as a 0/1 In this study the inverse protein folding problem
linear programming (01LP) problem, which is defined is reduced to the quadratical programming problem,
by a linear criterion to optimize over a set of Boolean which can be written as
variables under a conjunction of linear equalities and
inequalities. It is also represented as a problem of 0/1 minimize f0((cid:126)x)
(1)
quadratic programming (01QP), which has a quadratic subject to f((cid:126)x)≤0, i=1,2,...,N,
i
criteriontooptimizewiththesametypesofconstraints
as in 01LP. To solve these both models, the authors where
usedthecplexsolverandcomparedtheresultswiththe
CFN approach. fi(x)=(cid:126)xTAi(cid:126)x+2(cid:126)bTi(cid:126)x+ci,i=0,1,...,N.
In[6],theinversefoldingproblemissolvedusingge-
The initial problem is not convex, because the matri-
neticalgorithm. Thisalgorithmoptimizesthesequence
ces A are generally indefinite, therefore the problem is
of amino acids in a protein chain with a stochastic i
NP-hard and can not be solved efficiently. This paper
search similar to natural selection in evolution. This
studies semidefinite relaxation , Lagrangian relaxation
approach was later improved in [8] by expressing the
[3]andcontinuousrelaxationoftheproblem(1)tocon-
problemasahighlymulti-modaloptimizationproblem.
vex one, which can be solved in practical time. The
This is achieved by using a diversity measure as the
solution of the relaxed problem provides a lower bound
objective function through multi-objectivization. The
of the optimal value of (1). We use the rounding of
common feature of all genetic algorithms is that they
therelaxedsolutionto find anupper boundandan ap-
rely on evolutionary information of existing structures.
proximate the optimal value. This paper also suggests
Investigation [7] studies the Answer Set Program-
SequentialQuadraticProgramming(SQP)[4]andSim-
ming (ASP) approach for the inverse folding problem.
ulated Annealing approaches to solve (1). Finally, we
It uses DEE in combination with branch and bound al-
compare these relaxations and techniques and estimate
gorithmstoeliminatehighnumberofnon-optimalpairs
the efficiency and the quality of primary structure pre-
(amino acid, conformation), as in [1]. Then the prob-
diction of using different approaches on real proteins
lem is represented in the ASP model and solved with
structures from the SCWRL4 test set [5].
solver clingo4.
All the aforementioned investigations determined
2. Problem statement
both amino acid and conformation for all positions in
the polypeptide chain. These approaches used energy
Let the protein chain consist of N amino acids.
functions which consider interactions between all pos-
The set C ={1,2,...,20} contains indexes that encode
sible pairs of rotamers, for example, energy function,
all 20 possible amino acids. Let (cid:126)y = (y ,...,y ) be
implementedintheCPDdedicatedtoolosprey2.0[17]. 1 N
the sequence of residues, where y ∈C. Let functions
Thementionedworkscomparetheefficiencyofdifferent i
E : C2 →R define symmetrical pairwise interaction
algorithms, i.e. the time these algorithms need to find kl
energiesbetweenresiduesk andl. Thentheproblemof
the minimum of energy for given instances. The CFN
energy minimization is represented as
model demonstrated the best results overall, although
even it could not solve all the instances in reasonable
N N
time. However,thequalityofpredictionoftheprimary ∑∑E (y,y )→ min . (2)
ij i j
structure (the coincidence between the found structure j=1i=1 y1,y2,...,yN∈C
and the native one) is not investigated in these works.
The problem can also be written in the following way:
This quality depends on the choice of energy potential,
and there is a challenging problem to find a practical
minimize (cid:126)xTQ(cid:126)x
scoring potential, which does not always match with
(cid:126)x=[(cid:126)xT,...,(cid:126)xT]T
1 N
the true potential function [18]. subject to (cid:126)x ∈{0,1}20, k=1,...,N, (3)
The main distinction of our investigation is that k
(cid:107)(cid:126)x (cid:107) =1, k=1,...,N,
onlytheaminoacidsintheproteinchainarepredicted, k 0
and spatial conformations (rotamers) are not consid-
where
ered. This approach decrease the computational com-
plexity of the problem significantly. It allows to re-
 
[E ] [E ] ··· [E ]
ceive basic information about the primary structure of 11 12 1N
athpersoetqeuinenucseinogfmamucinholeascsidresseovuerncefso,ratnhdeamlloowstsctoomfipnlid- Q=[E...21] [E...22] .·.··. [E2...N],
catedinstances. Thus, thisstudyusesenergyfunctions [E ] [E ] ··· [E ]
N1 N2 NN
2
 
Eij(c1,c1) Eij(c1,c2) ··· Eij(c1,c20) Algorithm 1 Greedy algorithm
Eij(c2,c1) Eij(c2,c2) ··· Eij(c2,c20) Require:(cid:126)y =(y0,...,y0), y0∈C, k=1,N
Eij= ... ... ... ... , Ensure:(cid:126)yprsetdar=t (y1,1...,yN)N k
E (c ,c ) E (c ,c ) ··· E (c ,c ) repeat
ij 20 1 ij 20 2 ij 20 20
num changed←0
Q∈R20N×20N, E ∈R20×20.
ij for all k∈{1,...,N} do
yˆ :=argminf(y ,...,y ,y(cid:48),y ,...,y )
k 1 k−1 k k+1 N
Finally, it is practical to re-write the problem as y(cid:48)k∈C
if yˆ (cid:54)=y then
k k
minimize (cid:126)xTQ(cid:126)x y :=yˆ
k k
(cid:126)x∈{0,1}20N (4) num changed:=num changed+1
subject to A(cid:126)x=1N, until num changed>0
where
 
1···1 0···0 ······ 0···0 Oneachiteration,thisalgorithmsearchesforapro-
0···0 1···1 ······ 0···0 tein conformation of a lower energy than the current
A= ... ... ... ... , A∈{0,1}N×20N. one among conformations in the neighborhood of the
  current conformation. On each step it determines for
0···0 0···0 ······ 1···1
(cid:124)(cid:123)(cid:122)(cid:125) (cid:124)(cid:123)(cid:122)(cid:125) (cid:124)(cid:123)(cid:122)(cid:125) each position of the chain whether there exists such an
20 20 20 amino acid yˆ that would decrease the energy of the
k
Toestimatethequalityofprimarystructurepredic- protein by substitution y with yˆ . It stops when there
k k
tion, we use the scoring matrix BLOSUM62 [9]. This is no position to change, so this algorithm returns at
matrix is a substitution matrix used for calculating the least the local minimum of the problem (2).
degree of matching between two protein sequences.
Let B represents the BLOSUM62 matrix. Then, 3.2. Continuous relaxation
the score B(y,yˆ ) is the score for the substitution of
i j
residue y for yˆ at j-position in the chain. If the score ThematrixQin(4)couldberelaxedtosemidefinite
i j
is positive, the residues are interchangeable (or equal). positive by shifting its spectrum. Denote
Otherwise, this substitution is less likely to occur, and
the prediction is wrong. Let (cid:126)y ,(cid:126)y ∈ CN be na- Qˆ =Q−λmin(Q)·I20N, c=λmin(Q)·N. (6)
nat pred
tive and predicted sequences of length N, respectively.
Then
Then, the quality function can be defined as
(cid:126)xTQˆ(cid:126)x=(cid:126)xT(Q−λ I )(cid:126)x=(cid:126)xTQ(cid:126)x−λ (cid:126)xT(cid:126)x.
S((cid:126)y ,(cid:126)y )= ∑Nk=1B(((cid:126)ynat)k,((cid:126)ypred)k) (5) min 20N min
nat pred ∑Nk=1B(((cid:126)ynat)k,((cid:126)ynat)k) Considering(cid:126)x∈{0,1}20N andA(cid:126)x=1N,wehave(cid:126)xT(cid:126)x=N.
Therefore, optimization problem (4) can be written as:
AgoodqualitypredictioncorrespondstoS>0,andthe
best predictions have values of S close to 1. Negative minimize (cid:126)xTQˆ(cid:126)x+c
values of S mean poor quality predictions. (cid:126)x∈{0,1}20N (7)
The aim of this study is to compare different op- subject to A(cid:126)x=1 .
N
timization approaches to solve the problem (4). The
value of S is considered as the quality of the primary Now the matrix Qˆ is semidefinite positive due to (6),
structureprediction,whiletheoptimizationpowerofal- because all its eigenvalues are nonnegative. Consider
gorithms is estimated from the ability of the algorithm continuous relaxation of (7) by relaxing constraints(cid:126)x∈
topredictthesequencewiththelowestpossibleenergy. {0,1}20N to(cid:126)x≥(cid:126)0 :
20N
3. Algorithms and relaxations minimize (cid:126)xTQˆ(cid:126)x+c
(cid:126)x≥(cid:126)020N (8)
subject to A(cid:126)x=1 .
3.1. Greedy algorithm N
The problem (8) is convex now, and therefore can be
Given a problem (2) solved efficiently. The solution of (8) also gives the
lower bound for the optimum of the primary problem
N N
∑∑E (y,y )→ min , (4). An approximate solution of (4) and an approxi-
ij i j
j=1i=1 y1,y2,...,yN∈C mate optimum are found from(cid:126)x=[x(cid:126)1T,...,x(cid:126)NT]T as fol-
lows:
considerthegreedysearchalgorithmtofindanapprox- y =argmax xj, k=1,...,N.
imate solution: k k
j=1,20
3
Thesolutioncouldbealsoimprovedbyusingthegreedy The dual function is then
algorithm 1 with the solution of (8) as a starting posi-
g((cid:126)λ,(cid:126)u)= inf L((cid:126)x,(cid:126)λ,(cid:126)u)= inf (cid:126)xTP((cid:126)λ)(cid:126)x−(cid:126)qT((cid:126)λ,(cid:126)u)(cid:126)x−r((cid:126)u)=
tion. The final approximation will be the upper bound
(cid:126)x∈R20N (cid:126)x∈R20N
for the optimum value of the initial problem (4). −1(cid:126)qT((cid:126)λ,(cid:126)u)P+((cid:126)λ)(cid:126)q((cid:126)λ,(cid:126)u)−r((cid:126)u), if P((cid:126)λ)(cid:23)0 and
 4
3.3. Semidefinite relaxation (SDP) = (cid:126)q((cid:126)λ,(cid:126)u)⊥KerP((cid:126)λ),
−∞, otherwise.
SDP relaxation is described in [2]. Applying it to
the problem (4) we have Finally, the Lagrangian relaxation of (4) is:
maximize γ−r((cid:126)u)
minimize Tr(QX) (cid:126)λ∈R20N,(cid:126)u∈RN
X∈R20N×20N
(cid:126)x≥(cid:126)020N subject to γ ≤0, (10)
subject to X(cid:23)(cid:126)x(cid:126)xT, (cid:34) P((cid:126)λ) 1(cid:126)q((cid:126)λ,(cid:126)u)(cid:35)
2 ∈S20N+1.
Xii=xi, i=1,20N, 1(cid:126)qT((cid:126)λ,(cid:126)u) −γ +
2
A(cid:126)x=(cid:126)1 ,
N
3.5. Sequential Quadratic Programming
where Sn ={X∈Rn: X=XT(cid:23)0} is the set of sym-
+
metric semidefinite positive matrices of dimension N. Sequential Quadratic Programming (SQP) [4, 11]
By Schur’s complement [10], the problem is equivalent is one of the most effective methods for constrained
to nonlinear optimization. The idea of SQP is to solve
nonlinear problem using a sequence of quadratic pro-
minimize Tr(QX)
X∈R20N×20N gramming (QP) subproblems. The objective function
(cid:126)x≥(cid:126)020N of the subproblem on each step is a convex quadratic
(cid:20) (cid:21)
subject to X (cid:126)x ∈S20N+1, (9) approximation of the Lagrangian function.
(cid:126)xT 1 + First,werelaxtheconstraintsoftheinitialproblem
X =x, i=1,20N, (4)inordertoapplytheSQPmethod. Theconstraints
ii i
x ∈{0,1} should be substituted for 0≤x ≤1:
A(cid:126)x=(cid:126)1 . i i
N
minimize (cid:126)xTQ(cid:126)x
Again,therelaxedproblemisconvexandcanbesolved (cid:126)x∈R20N
efficiently using the methods of convex programming. subject to A(cid:126)x=1 , (11)
N
Theoptimumof (9)givesthelowerboundfor(4). The
0 ≤(cid:126)x≤1
20N 20N
binaryvector(cid:126)xandtheupperboundfortheinitialprob-
lem(4)arerecoveredfromthesolutioninthesameway The Lagrangian of (11) is
as in the previous case (continuous relaxation).
L((cid:126)x,(cid:126)λ,(cid:126)µ,(cid:126)u)=(cid:126)xTQ(cid:126)x−(cid:126)λT(cid:126)x+(cid:126)µT((cid:126)x−1 )+(cid:126)uT(A(cid:126)x−(cid:126)1 ).
20N N
3.4. Lagrangian relaxation (12)
Now let us formulate the QP subproblem for the k-th
TheLagrangianrelaxation[3,10]isagoodmethod step of the algorithm:
for getting a cheaply computable lower bound on the
1
optimal value of a non-convex quadratic optimization minimize d(cid:126)TH d(cid:126)+2(cid:126)x TQd(cid:126)
k k
problem. This method uses the fact that the dual of a d(cid:126)∈R20N 2
problem is always convex, and the weak duality guar- subject to A(cid:16)x (cid:126)+d(cid:17)=1 , (13)
k N
antees the optimum value of a dual problem to be the
lower bound to the initial one. 020N ≤xk(cid:126)+d≤120N.
The Lagrangian of (4) is
HereH isapositivedefiniteapproximationoftheHes-
k
sian matrix of the Lagrangian function (12), and (cid:126)x is
k
20N
L((cid:126)x,(cid:126)λ,(cid:126)u)=(cid:126)xTQ(cid:126)x+∑λ(x2−x)+(cid:126)uT(A(cid:126)x−(cid:126)1 )= a current iterate. The solution of (13) is used to form
i=1 i i i N a new iterate ((cid:126)xk+1,(cid:126)λk+1,(cid:126)µk+1,(cid:126)uk+1).
=(cid:126)xT(cid:16)Q+D((cid:126)λ)(cid:17)(cid:126)x−(cid:126)λT(cid:126)x+(cid:126)uT(A(cid:126)x−(cid:126)1 )= Each iteration of the algorithm consists of three steps:
N
    1. Updating the Hessian Matrix
=(cid:126)xTQ(cid:124) +(cid:123)D(cid:122)((cid:126)λ(cid:125))(cid:126)x−(cid:126)λ(cid:124)T−(cid:123)(cid:122)(cid:126)uTA(cid:125)(cid:126)x−(cid:126)u(cid:124)T(cid:123)(cid:126)1(cid:122)N(cid:125). 2. Solving a Quadratic Program
P((cid:126)λ) (cid:126)qT((cid:126)λ,(cid:126)u) r((cid:126)u) 3. Line Search and Merit Function
4
1. At each iteration a quasi-Newton approximation of states–allmannerofpossiblesequences. ThesetN(x)
the Hessian ∇2 L((cid:126)x ,(cid:126)λ ,(cid:126)µ ,(cid:126)u ) is calculated us- is a set of neighbors to the state x∈S. The energy
xx k+1 k+1 k+1 k+1
ingtheBroyden-Fletcher-Goldfarb-Shanno(BFGS)[12] function F(x) is the objective function. x0 ∈S is the
method. Denote starting position. Then the iteration of SA algorithm
is:
(cid:126)s =(cid:126)x −(cid:126)x ,
k k+1 k 1. Randomly choose the state from the neighbors of
(cid:126)y =∇ L((cid:126)x (cid:126)λ ,(cid:126)µ ,(cid:126)u )−∇ L((cid:126)x ,(cid:126)λ ,(cid:126)µ ,(cid:126)u ) a current state xk+1∈N(xk)
k x k+1 k+1 k+1 k+1 x k k k k
2. Assess the energy of the chosen state F =
Then the BFGS approximation for the next iteration k+1
F(x )
has the form k+1
Hk+1=Hk+(cid:126)(cid:126)yykT(cid:126)(cid:126)ysTk −H(cid:126)skT(cid:126)sHk(cid:126)sTk(cid:126)sHTk. (14) 3. cCeopmtpsaelreectFekd+1maonvde.FkIfaFnkd+1d<eciFdke,,twhehemthoevretios aacc--
k k k k k cepted. Else, the move is accepted with the prob-
(cid:18) (cid:19)
This approximation keeps Hk positive definite during ability exp −Fk+1−Fk .
the iterations. T(k)
2. With iterate ((cid:126)x ,(cid:126)λ ,(cid:126)µ ,(cid:126)u ) and a quasi-Newton ap- Inthisinvestigationproblem, thesetSisthesetof
k k k k
allpossiblesequencesofaminoacidsforagivenlength,
proximation of Hessian H , the subproblem QP (13) is
k
neighbors N(x) for the sequence x are the sequences,
formulated. ConsideringtheconditiononH tobepos-
k
in which 5 or less residues differ from x. The energy
itivedefinite,QPisconvexandcanbesolvedefficiently.
function F(x) is the energy of the protein structure.
New iterates for dual variables(cid:126)λ ,(cid:126)µ ,(cid:126)u are de-
k+1 k+1 k+1
terminedfromthesolutionimmediatelyasthevaluesof
4. Solution
thedualvariablesof (13)intheoptimalpoint. Thenew
iterate for the primal variable is set to(cid:126)x =(cid:126)x +α d(cid:126) ,
where d(cid:126) is the solution of QP (13), andk+t1he okptimkumk This section describes the process of solving the
k
problem (4) using methods, presented in the previous
length of the step α is determined at the next step.
k
section. First, we build the energy matrix for a pro-
3. Thesteplengthparameterα isdeterminedinorder tein structure. In the current study, we use two coarse-
k
toproduceasufficientdecreaseinameritfunction. We grained distance-dependent potentials. One of them is
use the following function: described in [14]. Another is DFIRE-Ca potential [15],
whose parameters are kindly provided by the authors.
(cid:16) (cid:17)
Ψ((cid:126)x)=(cid:126)xTQ(cid:126)x+(cid:126)wT(cid:126)x++(cid:126)pT 120N(cid:126)−x˜ +(cid:126)rT(A(cid:126)x−1N), They both require only the distance between two Cα
+ carbon atoms in a pair of amino acids to estimate the
(15)
energyofinteractionbetweentheseresidues. So,byex-
Here, if (cid:126)m=[(cid:126)wT,(cid:126)pT,(cid:126)rT]T, then (cid:126)m is calculated as
tracting coordinates of C atoms of the backbone of a
α
(cid:26) ((cid:126)m ) +π (cid:27) givenprotein, wecanbuildtheenergymatrix. Thatis,
k i i
(cid:126)mi=((cid:126)mk+1)i=max πi, , i=1,...,41N, for each pair of positions i,j in the chain, we calculate
2
the distance d between C atoms at these positions.
ij α
wherevector(cid:126)π=(cid:104)(cid:126)λT,(cid:126)µT,(cid:126)uT(cid:105)T isthevectorofdualvari- Then,theenergyoftheinteractionbetweenall210pairs
ofaminoacidslocatedatdistanced areestimated,and
ij
ables of QP (13). The parameter α is calculated then
k this is how matrices E are obtained. By doing so for
ij
from the minimization of the univariate function
all i,j∈1,N, we build the matrix Q.
Then different relaxations are applied to the prob-
Ψ((cid:126)x +α d(cid:126) )→min.
k k k αk lem (4):
• Lagrangian relaxation (10) is used to receive the
The problem (11) is then reduced to a sequential lower bound for the optimum value. It does not
solving of convex quadratical optimization problems. provide an approximate solution to the problem,
because there is no strong duality between (4) and
3.6. Simulated Annealing its dual.
• Greedy algorithm 1 is initialized with a random
Simulated Annealing (SA) [13] is used to compare
starting position.
different approaches. It is a probabilistic method, and
although it is unlikely to find the global optimum solu- • Continuous and SDP relaxations find both lower
tion, it may often find a solution very close to it. bound and approximate optimum of the problem.
Let T :Z →(0;+∞) be the exponentially decreas- We improve received solutions with the greedy al-
+
ingfunctionofthetemperature, andSbeafinitesetof gorithm in the next step.
5
• SQP and Simulated Annealing give an approxi-
mate optimum of the energy. 0
Finally,weestimatethequalityofoptimizationandpre-
diction for each method at each test protein structure. alue −5000
v
Thequalityofpredictionisobtainedfrom(5)usingthe m
u
BLOSUM62 matrix. m−10000
pti
o
e
5. Computational details mat−15000 Native energy
xi
o
pr Approximate optimum:
All the tested algorithms were implemented in Ap−20000 ContGreedy
RandGreedy
Python, using the CVXPY package [16] to solve dif- Anneal
SQP
ferentconvexproblemsthatareobtainedafterthepro-
posed relaxations. For the SQP algorithm, we used −250002pet1pxz6xia1atz1hm51w5r2b0j2ium1koe1wka1j231t1j1i4j2otu2pv22ibl1mkk1qkd1f941b2p1b9w1bgf1byi1c021f601fcq1fo91fpo1ftr1fvk1g611g8a1g8q1b8z1c481aho1bds
Proteins
its MATLAB implementation from the Global Opti-
mization Toolbox. The computation experiment was
Figure 1. Approximate energy optimum for different
set on a quad-core Intel Core(TM) i7-4700HQ CPU @
relaxations computed on the test set.
2.40GHz PC with 12 GB of RAM running on Ubuntu
14.04 LTS OS.
6. Results and Discussion 104
103
To compare the quality of different relaxation ap-
proaches in solving non-convex quadratic problems, we
e102
m
performedaseriesofcomputationalexperiments. More d ti
precisely, we used six algorithms, described in the sec- me101
tion 3, to find the minimum of the protein energy. The su
n
computational experiments were set on a test set of Co100 ContGreedy
protein structures extracted from the SCWRL4 bench- RandGreedy
Anneal
mark. The average length of sequences for these pro- 10-1 SQP
SDP
teins was 110 amino acids, so the dimension of(cid:126)x was
LagrRelax
about 2000. 10-20 20 40 60 80 100 120 140 160
Length of chain
Figure 1 presents the results of optimization the
energy using the potential from [14].
From this figure we can see that among the com- Figure 2. Computational time as a function of the
pared methods, SQP and SA demonstrated the best length of a protein chain for all the algorithms in
power of optimization. They returned very proximate logarithmic scale.
optimumenergyvaluesandpredictedsimilarsequences
of amino acids. The continuous relaxation approach
minimizes the objective function worse than SA and timization methods to solve the problem with cropped
SQP, but still returns the sequence with the energy, backbones, and the average time for each length was
lower than the native. We should note that the SDP measured. Also, we computed the mean of quality of
and Lagrangian relaxations methods cannot be applied predictionforsuchsequences. Thedependencesofcom-
to large optimization problems. This is because the putational time and quality on the length of a protein
number of their variables is quadratic with respect chain are shown in Figs. 2 and 3, respectively.
to the length of vector (cid:126)x, and also because the com- Figure 1 clearly demonstrates that all the relax-
putational requirements for these approaches exceed ations and all the tested methods found sequences of
the available resources for the sequences of a practi- amino acids with the energy lower, than the one of
cal length. For example, for N =100 the number of the native state. The greedy algorithm with a native
variables for these methods is more than 2·106. sequence of amino acids as the starting position also
For all the techniques and all the test structures, returned the sequence with a lower energy for all the
we measured the time spent on finding the solution of tested potentials. This means that the native struc-
the optimization problem. It was done in the following tures are not even the local minimum for the potential
way – for different lengths of sequences from 5 to 155, energyfunctionusedinthisexperiment. Therefore,the
we truncated the length of the structures of proteins twotestedpotentialsdonotsatisfythecrucialassump-
this upper limit. Then, we applied all the tested op- tion, i.e., that the native conformation has the lowest
6
structures that were found using the introduced tech-
−0.10 niques, had much lower energy than the native one,
ContGreedy
RandGreedy however the native and the predicted sequences coin-
Anneal
−0.15 SQP cided very weakly. We strongly believe that in order
SDP
to achieve better prediction results, an improvement of
y
alit−0.20 the potential functions for the inverse folding problem
qu is needed.
on −0.25
cti
di 8. Acknowledgments
e
Pr−0.30
This work was supported by RFBR, grant
−0.35
16-37-00111.
−0.400 20 40 60 80 100 120 140 160
Length of chain References
[1] David Allouche et. al. Computational protein design as
Figure 3. Dependence of prediction quality on the
an optimization problem. Artificial Intelligence, 2014.
length of a protein chain for all the tested methods.
[2] Stephen Boyd and Lieven Vandenberghe. Semidefinite
Programming Relaxations of Non-Convex Problems in
Control and Combinatorial Optimization. Communica-
energybetweenallpossiblestatesforagivenfixedback-
tions, Computation, Control, andSignalProcessing,pp.
bone, and thus the quality of primary structure predic-
279–287, 1997.
tion is extremely poor. In searching for a minimum of
[3] Alexandre d’Aspremont and Stephen Boyd. Relax-
energy,thealgorithmsselectresiduetypeswiththelow-
ations and Randomized Methods for Nonconvex QC-
est contact energies. This happens to be observed for
QPs. EE392o Class Notes, 2003.
pairs of cysteine-cysteine and tyrosine-tyrosine. As the [4] Philip E Gill, Walter Murray, and Michael A Saun-
result, the predicted structure often consist fully either ders. SNOPT:AnSQPAlgorithmforLarge-ScaleCon-
of cysteine residues, or of tyrosine residues. strainedOptimization. SIAM Journal on Optimization,
vol. 12, no. 4, pp. 979–1006, Jan. 2002.
7. Conclusion [5] Georgii G. Krivov, Maxim V. Shapovalov, and Roland
L. Dunbrack. Improved prediction of protein side-chain
conformations with scwrl4. Proteins: Structure, Func-
In this paper we proposed to formulate the Inverse
tion and Bioinformatics, 77(4):778–795, 2009.
Protein Folding Problem as a quadratic optimization
[6] David T. Jones. De novo protein design using pairwise
problem, where the objective function is the energy of potentialsandageneticalgorithm. ProteinScience,vol.
the protein, which is assumed to have the lowest value 3, no. 4, pp. 567–574, Dec. 2008.
forthenativeconformation. Theresultsshowsthatthe [7] Jacques Nicolas Hugo Bazille. Computational protein
SQP and SA approaches have the best optimization design: trying an answer set programming approach to
power on the test set of protein structures extracted solve the problem. 10th Workshop on Constraint-Based
from the SCWRL data set. These two methods Methods for Bioinformatics, 2014.
[8] Wiktor Jurkowski Sune S. Nielsen, Gregoire Danoy.
returnedverysimilarapproximateoptimumsfornearly
A Novel Multi-objectivisation Approach for Optimising
all the structures. The SDP approach turned out to
theProteinInverseFoldingProblem.EvolutionaryCom-
increase the computational cost of the relaxed problem
putation, Machine Learning and Data Mining in Bioin-
rapidlyandrequiredunreasonabletimeforthesolution
formatics, 2015.
for any practical case. Continuous relaxation gave the
[9] Henikoff JG. Henikoff S. Amino acid substitution ma-
worse approximation for all the test cases, compared trices from protein blocks. Proceedings of the National
to SQP and SA, but consumed less time to compute Academy of Sciences, vol. 89, no. 22, pp. 10915–10919,
the answer. All materials needed for conducting an Nov. 1992.
experiment can be found at https://sourceforge. [10] Lieven Vandenberghe Stephen Boyd. Convex Opti-
net/p/mlalgorithms/code/HEAD/tree/Group374/ mization. Cambridge University Press, 2004.
Ryazanov2016InverseFolding/code/. [11] J. Nocedal and S. J. Wright. Numerical Optimization,
SecondEdition. SpringerSeriesinOperationsResearch,
Although the tested algorithms demonstrated a
2006. Chapter 18.
good capacity in solving the optimization problem, the
[12] R. Fletcher. A new approach to variable metric al-
actual quality of the primary structure prediction was
gorithms. The Computer Journal, vol. 13, no. 3, pp.
poor. This is mainly because the potential functions
317–322, Mar. 1970.
that we used for the energy estimation do not follow [13] Dimitris Bertsimas, John Tsitsiklis, et. al. Simulated
theassumptionthatthenativesequencehasthelowest annealing. Statistical science, vol. 8, no. 1, pp. 10–15,
energy for a given geometrical shape. Therefore, the Feb. 1993.
7
[14] C. A. Floudas R. Rajgaria, S. R. McAllister. A novel
high resolution C −C distance dependent force field
α α
based on a high quality decoy set. Proteins: Struc-
ture, Function, and Bioinformatics, vol. 65, no. 3, pp.
726–741, Nov. 2006.
[15] Liu S et al. Zhang C. An accurate, residue-level, pair
potential of mean force for folding and binding based
onthedistance-scaled,ideal-gasreferencestate. Protein
Science, vol. 13, no. 2, pp. 400–411, Feb. 2004.
[16] Steven Diamond and Stephen Boyd. CVXPY: A
Python-Embedded Modeling Language for Convex Op-
timization.JournalofMachineLearningResearch,2016.
[17] P. Gainza, K.E. Roberts, I. Georgiev, R.H. Lilien,
D.A. Keedy, C.Y. Chen, F. Reza, A.C. Anderson, D.C.
Richardson, J.S. Richardson, et al. Osprey: protein de-
signwithensembles,flexibility,andprovablealgorithms.
Methods Enzymol, 2013.
[18] Ting Lan Chiu and Richard A.Goldstein. Optimizing
potentials for the inverse protein folding problem. Pro-
teinEngineeringDesignandSelection,vol.11,no.9,pp.
749–752, Sep. 1998.
8