Transcript Protein design as an inverse problem

Computational Protein Design:
A problem in combinatorial
optimization
CSE 549 Guest Lecture
September 17, 2009
David Green
Applied Mathematics & Statistics
What is a protein?
• Polymers (chains) of amino acids.
– There are 20 different amino acids that can be
part of the chain.
• Machines of the cell.
– It’s proteins that do most of the work involved in
life!
Polymers of amino acids.
• Amino acids link to form polypeptides.
– There is a backbone of constant composition.
– There are side chains that vary.
The twenty amino acids.
• AA side chains vary
from:
– Big to small.
– Non-polar (all C and
H) to polar.
– Positive to
negative.
– Flexible to rigid.
The machinery of life.
• Protein sensors (receptors) are responsible for all the
senses (sight, smell, taste, touch, hearing).
• Enzymes are proteins the catalyze chemical
reactions, like the ones that convert food to energy.
• Specialized structural proteins make skin elastic, and
make the lens of the eye work.
• Muscles are primarily composed of proteins that
combine structural and enzymatic parts to make a
machine.
Why design proteins?
• New sensors based on biology.
– Proteins have been engineered to detect TNT (explosive)
and sarin (nerve gas).
• Proteins are used as treatments for many diseases.
– Protein engineering has helped improve proteins that are
given to cancer patients on radiation or chemo-therapy.
– Work in the Green lab is on-going to design proteins for
use as anti-HIV prophylatics.
• Many nanotechnology applications that haven’t even
been considered yet!
Where do proteins come from?
• The genome contains instructions for every protein
in a cell.
– A few HUGE molecules of DNA.
• Each gene is the code for one protein.
– There are ~30,000 genes in humans.
• Genes are expressed through an intermediate
molecule, RNA.
– Many copies of each protein can be made.
The Central Dogma of Molecular
Biology.
• Then proteins do the work!
How do proteins work?
• Proteins fold into a unique 3-dimensional
structure.
• The amino acid sequence of a protein dictates
it’s structure.
• The function of a protein is controlled by it’s
structure.
Many polymers are long, unstructured
chains.
• Polyethylene
– Is made of long
chains of the same
monomers.
– Adopts a random
mesh of interweaving strands.
– This structure gives
us PLASTIC!
DNA has the same structure for every
sequence.
• The “double-helix” is a great structure for storing
and replication information.
Protein structures are well-defined
and diverse!
• One chain or
many.
• Elongated or
globular.
• Many forms of
symmetry (or
none).
What does a protein look like?
• Cyanovirin – A protein that inhibits the
entry of HIV into human cell.
What does a protein look like?
• The atoms of a protein form a compact,
well-packed cluster.
What does a protein look like?
• A protein can be thought of as a nearly
solid object.
What does a protein look like?
• Simplified cartoons make the structure
easier to see.
What does a protein look like?
• The path of the backbone of a protein is
called it’s “fold”.
What does a protein look like?
• Different types of amino acids are found all
along the protein chain.
What does a protein look like?
• Each amino acid has a side chain that
protrudes from the backbone.
What does a protein look like?
• Many proteins bind other molecules, like
the sugar molecules here.
What does a protein look like?
• Binding interfaces are usually a close fit of
two complementary surfaces.
What does a protein look like?
• The core of a protein is key in keeping a
stable structure.
Many side chains fill the core.
The core is well packed …
… with groups from all along the chain.
Each side chain fits perfectly.
What is a protein?
• A protein is a complicated three-dimensional
structure, made up by an amazing 3-D jigsaw
puzzle of interlocking amino acids.
• Amino acids pack together not just
geometrically, but with complementary
chemical groups as well.
• Proteins move too, but we’ll ignore that for
now.
How can we design one?!?
1. Choose a fold (path of the backbone).
2. Pack the core with the right set of amino
acids to achieve the desired fold.
3. Choose other amino acids to achieve the
desired function (such as binding to a target
molecule, or getting the right molecular
motions).
Structure prediction is a forward
problem.
•Given a protein
sequence, what is the
structure that it will
adopt (fold to)?
•This is a VERY hard
problem, and it not yet
fully solved.
•Prediction is difficult
because you are stuck
with what nature gives
you.
Protein design is an inverse problem.
•Given a desired 3dimensional protein
structure, what is a
sequence that will fold
to that structure?
•We have the freedom
to add constraints that
simplify the problem.
•As a result, methods
for protein design have
had many successes.
• Pabo. Nature 301: 200 (1981).
• Drexler. PNAS 78: 5275-5278 (1981).
A designed sequence should fold
according to design.
•ANY sequence which
folds to the correct
target structure (and
carries out the desired
function) can be
considered a successful
design
•There is more than
one right answer, unlike
in prediction!
Choosing a backbone fold.
• The structure dictates the function, and a big
part of structure is the fold.
• We still don’t really know how to choose the
“best” fold.
• Instead, we just borrow from nature –
redesign a natural protein to do something
new.
Zinc finger proteins bind DNA.
A Zinc ion holds them together.
• The protein will not
fold if zinc is not
present.
• The protein only binds
DNA when it is folded.
• A group at Caltech set
out to design a zinc
finger that doesn’t
need zinc!
1997: The first fully automated protein
design!
• Dahiyat and Mayo. Science 287: 82-87 (1997).
Designing function.
• Making a molecule bind is like designing a the
core – we want to make the interface between
the two pieces complementary.
• Other functions are a lot trickier … and we
don’t have good ways to solve them yet, but
we’re on our way.
2003: A Duke group designs a set of
protein sensors.
• Looger, Dwyer, Smith and Hellinga. Nature 423: 185-190 (2003).
Protein design is a BIG problem.
• The zinc finger is one of the smallest protein domains …
about 30 amino acids long.
• How many different 30 amino acid polypeptides are
there?
– Choose from any of 20 amino acids at each position.
– Total sequences = 2030 = 1x1039
• Mass of earth = 6x1027 g
• Mass of a grain of sand ~ 1x10-3 g
– A billion earths’ worth of sand grains
• Enumeration of possible states is beyond impossible
— must take advantage of need to achieve
complementary interactions between amino acids.
Many different structures are possible.
• An arginine and a glutamate interact.
Many different structures are possible.
• An arginine and a glutamate interact.
Many different structures are possible.
• An arginine and a glutamate interact.
Many different structures are possible.
• An arginine and a glutamate interact in several different
conformations.
Really Big!!!
• Amino-acid side chains are flexible.
• But not every shape (conformation) is equal.
• Each amino acid has a set of preferred conformations
(rotamers).
– 1 to 80 per amino acid.
• Instead of choosing from 20 amino acids … we need
to choose from ~400 (at least) amino acid rotamers!
– Total structures = 40030 = 1x1078
• (approx. number of atoms in the universe!!!!!)
Packing side chains – a puzzle.
• How do you solve a jigsaw puzzle?
• Impossible to try all combinations of piece
placement
– Unique ways of placing N pieces on a grid is (4N)(N!)
• For N=100, (1.6x1060)(9.3x10157)= 1.5x10218
• Trying each piece one by one is better, but still
infeasible
– Number of iterative tries for a N piece puzzle is:
N
N
4
4( N  i )( 2i  2)  N ( N  4)( N  1) 

3
i 1
• For N=100, 1.37x106
 4( N  i)(2i  2)
i 1
N
N


 8 N 2  ( N  1) i   i 2 
i 1
i 1


N ( N  1) N ( N  1)( 2 N  1) 

 8 N 2  ( N  1)


2
6

4
  N ( N  4)( N  1) 
3
Packing side chains – be smart.
• How do you solve a jigsaw puzzle?
• Group pieces by colors and patterns.
• Iterate over matching of pieces that are
complementary
– Shape is important.
– The pattern must also match.
Pattern matching in proteins?
• What does it mean for two amino acids in the
core of a protein to “match”?
– Must fit close together (but not too close) 
Steric complementarity.
– Neighboring atoms must have complementary
charges (neutral likes neutral, positive likes
negative)  Electrostatic complementarity.
Steric fit: Lennard-Jones potential.
 r0 12  r0  6 
U  U o    2  
 r  
 r 
• Van der Waals attraction
between atoms at moderate
distances.
• Repulsion of atoms from
one another at short
distances.
• If atoms are not nearby, the
energy between them will
be very close to zero.
• The total score of the
“goodness” of fit in a
molecule is the sum of
the energy for every pair
of atoms.
Electrostatic fit: Coulomb’s Law
q1q2
U  kC
r
• Atoms in molecules can be thought
of as having tiny charges on them,
even if the total charge on a
molecule is zero.
• Coulomb’s Law describes the
energy of how two charges interact.
• The overall electrostatic fit is
calculated by adding up the energy
of all pairs of atoms.
– Like charges give a positive value.
– Opposite charges give a negative.
– Neutral (zero charge) groups don’t
matter.
The total energy describes the fitness
of a structure.
• Van der Waals + Coulomb’s Law, for every pair
of atoms, and all added up.
• Negative energies are favorable, positive
energies unfavorable.
– Nature works to MINIMIZE energy.
Position 1
Protein Design as a Discrete
Conformational Search
Position 2
Position 3
Conformational states of system
Tree Pruning with Dead End Elimination
Molecular Mechanics energy
 Van der Waals
 Coulombic
 Bond, angle, dihedral
i , j are positions in sequence
Ri is rotamer choice at
position i

E ( R) 
 Ei,Ri
i positions
1 E pair
Ei , Ri  Eiself

 2 (i,Ri ),( j ,R j )
, Ri
i j
The Dead-end Elimination Theorem
Given two rotamer choices X, and Y at position I, if the best energy
of X (with any choice of rotamers at other positions) is worse than
the worst energy of Y, then X can not be part of the global energy
minimum.
Need to make the comparison,
 ?

min
 ( E ( X i , R))  max
 ( E (Yi , R))
R
R
but the min and max functions require evaluating all states.
Making DEE feasible.

 self

pair
1
E ( X i , R)    E p , R p   2 E( p , R p ),( q , Rq )  * Ri is used to replace Xi in the first equation.
p 
q p

 self

pair
1 E pair
 Eiself

E

E

 p , R p  2 ( p , R p ),( q , Rq ) 


,Xi
( i , X i ), ( q , Rq )
q i
p i 
q  p ,i

Last sum is invariant with respect to our choice at position i, and thus:




 ?

 self
?
 self

pair
pair
min
E ( X i , R)  max
E (Yi , R)  min

  Ei , X i   E( i , X i ),( q , Rq )   max

  Ei ,Yi   E( i ,Yi ),( q , Rq ) 
R
R
R
R
q i
q i




But note that:
 self

pair
self
pair
min
  Ei , X i   E( i , X i ),( q , Rq )   Ei , X i   min E( i , X i ),( q , Rq )
Rq
R
q i
q i


 self

pair
self
pair
max
  Ei ,Yi   E( i ,Yi ),( q , Rq )   Ei ,Yi   max E( i ,Yi ),( q , Rq )
Rq
R
q i
q i






This gives a sufficient condition for the DEE theorem to hold:
self
i, X i
E

  min E
q i
Rq
pair
( i , X i ), ( q , Rq )
 E
?
self
i ,Yi

  max E(pair
i ,Yi ), ( q , Rq )
q i
Rq

min and max are evaluated over rotamers at a single position, so entirely feasible!
Improving the bounds for DEE
Original problem statement was to compare the best structure with Xi at
position I to the worst with Yi; it is a easier to satisfy, but still sufficient,
criterion to find the single set of choices at all other positions with the minimum
difference between choice Xi and Yi.






 ?


 ?
min
E ( X i , R)  max
E (Yi , R)  min
E ( X i , R)  E (Yi , R)  0



R
E
self
i, X i
R
E
self
i ,Yi

R

?
pair
pair
 min
  E( i , X i ),( q , Rq )  E( i ,Yi ),( q , Rq )   0
R
 q i

Again, we use the same trick to bound the minimum in a feasible manner:



pair
pair
pair
pair
min
  E( i , X i ),( q , Rq )  E( i ,Yi ),( q , Rq )    min E( i , X i ),( q , Rq )  E( i ,Yi ),( q , Rq )
R
R
 q i
 q i q

This gives a alternate sufficient condition for the DEE theorem to hold:
E
self
i, X i
E
self
i ,Yi
   min E
q i
Rq
pair
( i , X i ),( q , Rq )
E
pair
( i ,Yi ),( q , Rq )
 0
?
This is a tighter bound on the “true” desired comparison, since:
 min E
q i
Rq
pair
( i , X i ), ( q , Rq )
  max E
q i
Rq
pair
( i ,Yi ),( q , Rq )
  min E
q i
Rq
pair
( i , X i ), ( q , Rq )
 E(pair
i ,Yi ), ( q , Rq )

DEE as an iterative algorithm
As rotamers are flagged as incompatible with the global minimum solution, the
min/max functions are evaluated over a smaller and smaller set of choices, and
so additional iterations of the comparison can eliminate more possibilities.
Thus, the algorithm can be outlined as:
While (any rotamers eliminated ) {
For each position in the sequence (i){
For each rotamer choice (X) at position i {
For each rotamer choice (Y ne X) at position i {
If DEE criterion is satisfied {
Eliminate choice X at position I
}
The order of each cycle is ~N*R2, where N is
}
the number of positions, and R is the number
}
of choices at each position. However, as R
}
decreases with iteration, each subsequent
}
cycle costs less.
DEE identifies branches that are
incompatible
with
the
global
minimum
Position 1
Position 2
Position 3
Conformational states of system
The pruned tree can be much smaller
Position 1
Position 2
Position 3
Conformational states of system
DEE can be highly effective
Example of a 5-position design, with 306 choices per position, done in a simple
MATLAB script:
Iteration
Iteration
Iteration
Iteration
Iteration
Iteration
Iteration
Iteration
Iteration
Iteration
Iteration
Iteration
Iteration
Iteration
Iteration
Iteration
Iteration
0:
0:
1:
1:
1:
2:
2:
2:
3:
3:
3:
4:
4:
4:
5:
5:
5:
306 306 306 306 306
Structures = 2.682916e+12
143 198 145 42 33
Structures = 5.690265e+09 from 2.682916e+12
Elapsed time is 96.461053 seconds.
52 83 55 11 8
Structures = 20889440 from 5.690265e+09
Elapsed time is 11.604577 seconds.
4 12 36 6 5
Structures = 51840 from 20889440
Elapsed time is 0.625266 seconds.
4 12 35 6 5
Structures = 50400 from 51840
Elapsed time is 0.148598 seconds.
4 12 35 6 5
Structures = 50400 from 50400
Elapsed time is 0.045631 seconds.
DEE: Notes and Caveats
Dead-end elimination is not guaranteed not to eliminate any
choices … in this case computational expense is used at zero gain.
However, experience suggests that in the case of protein design,
the algorithm is highly efficient.
For large design problems, even a highly efficient pruning can leave
a tree which is too large to be searched by enumeration (such as
depth-first search); for example, consider an original space of 10100
states, reduced to 1020.
An efficient bounding heuristic can be defined using similar tricks as
discussed here; this can be used in the A* algorithm to find the
global minimum within the remaining space.
The DEE criterion can be extended to pairs, triples, and n-tuples of
positions. Most applications use only singles and pairs.
Optimization as a tree search.
How to choose a path through a tree, which the goal of
reaching the global minimum state?
First, define an energy to be associated with each step through
the tree:
i 1
pair
Ei , Ri  Eiself

E
 (i,Ri ),( j ,R j )
, Ri
j 1
This is the energy of placing rotamer choice R at position i, given that
all positions in the tree above i have been selected.
Note this is similar, but not identical, to the definition used with DEE.
At the leaf node, the state energy is then the sum of all individual
path energies:

E ( R) 
E
i , Ri
i positions
Thus, we wish to find the path with the lowest total.
Efficient search using a heuristic (A*)
Challenge is that we do not know the total until we
have traversed the tree to a leaf!
Assume we have some heuristic, that estimates the best possible
energy that we can achieve down some path. For choice R at
*
position I, we will call this E i ,R
i
When choosing a path to take at a given node, we know the path we
have already taken, and thus the true cost. We thus combine this
information with the heuristic in deciding what step to take:
i 1
f (i, Ri )  E *i ,Ri   E j , R j
j i
At a leaf node, the heuristic is 0, and this gives the total energy.
However, if we only use the heuristic, how do we know we get the
correct solution?
Guaranteed optimality with A*
1. Allow backtracking.
At every step, not only choose between the possible paths from the
current node, but also the paths from all nodes which have been
visited in the past.
In this way, if the heuristic turned out to be poor for a given path (and
the true energy became large), a new path is chosen.
2. Use a heuristic that bounds the true solution.
Consider a heuristic that is guaranteed to be lower than the true best
energy down a path (an optimistic prediction of the best energy). When
a leaf is reached, a comparison between the true energy of that leaf
and the heuristic energy of the un-followed choice can be made. Since
the heuristic is an optimistic guess, if the true energy is lower than the
heuristic for all other choices, it must be the global minimum.
Thus, it is possible to have a guarantee that the solution found is
the true global minimum!
Defining the bounding heuristic
Recall our energy definitions:
Ei , Ri  E

E ( R) 
i 1
self
i , Ri
E
j 1
pair
( i , Ri ),( j , R j )
The optimal heuristic is:
i 1
f opt (i, Ri )   E j , R j  min

j 1
R
E
i , Ri
i positions
N
E
j i
j,R j
Now, we bound the second term:
 self j 1 pair
 N
 self j 1 pair

min
   E j , R j   E( j , R j ),( k , Rk )    min
  E j , R j   E( j , R j ),( k , Rk ) 
R
R
j i 
k 1
k 1
 j i


N
 self
 self j 1 pair
 N 
 j 1 pair
 


min
E
  E j , R j   E( j , R j ), ( k , Rk )    min  E j , R j  min
 


(
j
,
R
),
(
k
,
R
)


j
k


R
R
j i
k 1

 j i  R j 
 k 1
 
N

 self
 j 1 pair
  N 
 self j 1
pair





min
E

min
E

min
E

min
E

 R  j,R j



( j , R j ),( k , Rk )    
j,R j
( j , R j ),( k , Rk )
Rj 
Rk
R 
j
j i 
k 1
 k 1
  j i 



N





Overview of A* for protein design
Initialize a list of traveled nodes with the root.
While (no mininum leaf in list) {
Select minimum f* of all paths from nodes in list.
Add this new node to list.
If (new node is leaf) {
Compare leaf energy to minimum f* from list.
If (Leaf energy < min(f*)) {
Leaf is global minimum.
}
}
Our final heuristic is:

j 1
j 1
N 



pair
self
pair

f * (i, Ri )    E self

E

min
E

min
E




j,R j
( j , R j ), ( k , Rk ) 
j,R j
( j , R j ),( k , Rk )
Rj 
Rk
j 1 
k 1
k 1
 j i 

i 1
The second term involves a minimum over each choice at
following position, which is order (N-i)*R, with an inner minimum
over the same, order (N-j)*R. Thus, the cost is less than (N-j)2*R2.




A*: Notes and Caveats
Performance of A* is not guaranteed — in the worst case, the
entire tree must be enumerated to find the solution. Again,
however, experience suggests that the algorithm is highly efficient
for protein design.
In many cases, we wish not only to have the global minimum, but
all solutions within a cutoff of the minimum.
A* can be adapted to solve this problem as well, but care must be
taken in the first step of pruning with DEE — the elimination
criterion must be modified to prevent elimination of low, but not
minimum, energy states.
References
1.
2.
3.
4.
5.
6.
7.
8.
Protein design as an inverse problem:
Pabo. Nature 301: 200 (1981).
Drexler. PNAS 78: 5275-5278 (1981).
Examples of successful protein design:
Dahiyat and Mayo. Science 287: 82-87 (1997).
Looger, Dwyer, Smith and Hellinga. Nature 423: 185-190 (2003).
Development of the Dead-End Elimination algorithm:
Desmet, DeMaeyer, Hazes and Lasters, Nature 356: 539-542 (1992).
Goldstein, Biophys. J. 66: 1335-1340 (1994).
Gordon and Mayo. J. Comput. Chem. 19: 1505-1514 (1998).
Formulation of A* for protein design:
Leach and Lemon, Proteins 33: 227-239 (1998).