Abstract

The success of antibody-based pharmaceuticals has led to a resurgence in interest in computational structure-based design. Most efforts to date have been on the redesign of existing interfaces. These efforts have mostly neglected the inherent flexibility of the receptor that is critical for binding. In this work, we extend on a previous study to perform a series of designs of protein binding interfaces by incorporating receptor flexibility using an ensemble of conformers collected from explicit-solvent molecular dynamics (MD) simulations. All designer complexes are subjected to 30 ns of MD in explicit solvent to assess for stability for a total of 480 ns of dynamics. This is followed by end-point free energy calculations whereby intermolecular potential energy, polar and non-polar solvation energy and entropy of ligand and receptor are subtracted from that of the complex and averaged over 320 snapshots collected from each of the 30 ns MD simulations. Our initial effort consisted of redesigning the interface of three well-studied complexes, namely barnase–barstar, lysozyme–antibody D1.3 and trypsin–BPTI. The design was performed with flexible backbone approach. MD simulations revealed that all three complexes remained stable. Interestingly, the redesigned trypsin–BPTI complex was significantly more favorable than the native complex. This was attributed to the favorable electrostatics and entropy that complemented the already favorable non-polar component. Another aspect of this work consisted of grafting the surface of three proteins, namely tenascin, CheY and MBP1 to bind to barnase, trypsin and lysozyme. The process was initially performed using fixed backbone, and more than 300 ns of the explicit-solvent MD simulation revealed some of the complexes to dissociate over the course of the trajectories, whereas others remained stable. Free energy calculations confirmed that the non-polar component of the free energy as computed by summing the van der Waals energy and the non-polar solvation energy was a strong predictor of stability. Four complexes (two stable and two unstable) were selected, and redesigned using multiple conformers collected from the MD simulation. The resulting designer systems were then immersed in explicit solvent and 30 ns of MD was carried out on each. Interestingly, those complexes that were initially stable remained stable, whereas one of the unstable complexes became stable following redesign with flexible backbone. Free energy calculations showed significant improvements in the affinity for most complexes, revealing that the use of multiple conformers in protein design may significantly enhance such efforts.

Introduction

The design of new proteins to perform new functions is of significant interest to the biomedical research community (Binz et al., 2005; Butterfoss and Kuhlman, 2006; Rosenberg and Goldblum, 2006; Lippow and Tidor, 2007; Razeghifard et al., 2007). Designer proteins can serve as antagonists or agonists to modulate cell signaling pathways. Agonists promote signaling and therefore could be valuable in various biomedical applications such as regenerative medicine. Antagonists that block signaling typically by blocking protein–protein interactions could be useful in the development of biopharmaceuticals. The success of therapeutics such as bevacezumab and etanercept has led to a resurgence in antibody-based biopharmaceuticals (Jain et al., 2007). However, antibodies have significant limitations, such as exorbitant manufacturing costs, intellectual property impediments, difficulty in accessing large tumors due to the size of the protein and difficulty in designing an antibody to act as agonists.

The increased interest in protein-based biopharmaceuticals has also spurred interest in structure-based computational design, which has been made more attractive in light of the significant amount of structural information that is now available, as well as increasingly sophisticated computational tools and algorithms (Binz et al., 2005; Rosenberg and Goldblum, 2006; Lippow and Tidor, 2007; Razeghifard et al., 2007). Although traditionally computational design has been focused on de novo design of stably folded proteins, an area of intense interest is the design of a new function through the modulation of the binding interface of a protein. A number of successes in this area have been reported in the literature (Vita et al., 1999; Zondlo and Schepartz, 1999; Montclare and Schepartz, 2003; Sia and Kim, 2003; Liu et al., 2004; Yin et al., 2007) such as the transfer of a binding epitope of an HIV-1 C-peptide onto the surface of a GCN4 leucine zipper as stable coiled coil (Sia and Kim, 2003), the enhancement of affinity of the antibody bevacezumab (Lippow et al., 2007), the creation of new DNase-inhibitor protein pairs (Kortemme et al., 2004) and finally the design of a mixed 21-residue miniprotein with mixed α/β structure (Ali et al., 2005). It is worth noting that one common feature of these studies is that the grafting involved was done for the purpose of strengthening its interaction with the receptor. Typical protein–protein interfaces, however, are large (between 1200 and 2000 Å2) such that multiple secondary structures within the binding partners participate in the binding (Keskin et al., 2008). These interactions present significantly greater challenges for the structure-based computational design of protein scaffold.

Although several successes in protein grafting have been reported (Marvin and Hellinga, 2001; Yang et al., 2005; Razeghifard et al., 2007; Jiang et al., 2008), there is only one example to our knowledge that involves the design of a protein surface to bind to another unrelated protein (Liu et al., 2007). The challenges arise from a number of factors including limitations in the performance of scoring functions, lack of explicit treatment of water molecules and the neglect of entropic effects. In particular, the lack of receptor flexibility during the design process has contributed to the difficulty in designing stable complexes due to a significantly narrower conformational space that is sampled with a single structure (Bonvin, 2006; Lensink and Mendez, 2008). In protein design, several studies have attempted to include backbone motion (Fu et al., 2007; Humphris and Kortemme, 2008). Kortemme and coworkers developed a design strategy such that the sequences are selected to conform to an ensemble of conformers using their ‘backrub’ approach; their predictions qualitatively match experimentally observed sequence preferences (Humphris and Kortemme, 2008). Keating and coworkers employed normal mode analyses to sample conformers that were used to design diverse helical BH3 ligands that bind to the anti-apoptotic protein Bclx (Fu et al., 2007). Kuhlman and Baker used backbone flexibility for the de novo design of Top7, a 93-residue α/β protein with a novel sequence and topology (Kuhlman et al., 2003). A number of studies have worked to incorporate receptor flexibility in protein docking as well (Smith et al., 2005; Bonvin, 2006; Krol et al., 2007; Chaudhury and Gray, 2008; Lensink and Mendez, 2008).

In this work, we extend on a previous study (Liang et al., 2009) by explicitly incorporating receptor flexibility into the design process and validating the resulting complexes through extensive explicit-solvent molecular dynamics (MD) simulations. Here, we expand on our previous work in several respects. First, we assess the performance of our design scoring function through participation in a CAPRI competition. Second, we test the scoring function for its ability to predict free energy change upon alanine mutation at the binding interface. Third, we employ an ensemble of conformers collected from MD simulations of the native complex to explicitly incorporate backbone motion in our interface design efforts. Finally, we explore the stability of all complexes that emerge from this effort through a series of extensive explicit-solvent MD simulations and free energy calculations that consisted of assessing the contributions from electrostatics, non-bonded energy, polar and non-polar solvation free energies and entropy calculations.

Methods

Protein docking

ZDOCK2.3 is a fast rigid-body protein–protein docking algorithm that was used to dock the scaffold protein onto the receptor protein (Chen et al., 2003). The scoring function used for docking combines pairwise shape complementarity with desolvation and electrostatics. It includes a filter feature to specify the residues that are involved in binding. Initially, the native complex structure was superimposed to the unbound structure of the receptor. Those residues of the unbound structure of receptor, which were 10 Å far away from the superimposed ligand, were assigned an unfavorable desolvation parameter of 19, as defined by ZDOCK. Thus, the scaffold proteins would be docked onto the binding site of the receptor protein even in global search. The default docking parameters were used. The top 10 conformations were selected for visual analysis. The structure with no backbone clashing and the smallest cavities between the binding partners was selected.

Interface design

The protein interface was defined as the set of surface residues with side chains possessing a solvent-accessible surface area that decreased upon complexation. The interface residues of ligand protein in the complex were mutated to form favorable interactions. Side-chain conformation of interface residues of the receptor protein were modeled but not mutated. Side chains of other residues and backbone conformations were fixed for both the receptor and ligand proteins. The following equation was used to calculate the interaction energy between the positioned side-chain rotamer, the typical side-chain conformation and other parts of the protein:  

(1)
formula
where SAcontact, Voverlap, EH-bond, Eelec, ΔSApho and ΔSAphi are atom contact surface, overlap volume, hydrogen-bonding energy, electrostatic interaction energy, buried hydrophobic solvent-accessible surface and buried hydrophilic solvent-accessible surface between the rotamer and other parts of the protein, respectively; Fphi the fraction of buried surface of non-hydrogen-bonded hydrophilic atoms; Vexclusion the normalized solvent exclusion volume around charged atoms; f1 the observed frequency of the rotamer; f2 the observed frequency of the amino acid given a backbone conformation; Nss-bond the flag of disulfide bridge (1 or 0); ΔGref the reference value of each amino acid. The weights of these energy items together with the reference values were optimized so that the native residue was predicted energetically favorable at each position of the training proteins.

The binding affinity was computed using the following equation:  

(2)
formula
where SAcontact(I) is the contact surface area between two proteins; ΔSApho(I) and ΔSAphi(I) the accessible surface area difference between the complex and the unbound binding partners. Both side- and main-chain atoms are included in the calculation. A reference value of the protein–protein interaction is added to Eq. (2) instead of amino acid reference value in Eq. (1) to calculate binding free energy. We hypothesized that the interface reference value is proportional to the buried solvent-accessible surface area of the interface (ΔSAburied). A strong correlation between E(I), ΔGexp(I), ΔSAburied was found by multiple linear regression analysis of native complexes and binding data:  
(3)
formula

The last term in Eq. (3) is interface reference value (Liang et al., 2007). We used Eq. (3) to select near-native conformations in ZDOCK2.3 decoy sets (Chen et al., 2003). In 20 out of 48 targets, near-native decoys were predicted as the lowest binding free energy. We also participate CAPRI (Critical Assessment of PRediction of Interaction) (Janin et al., 2003) using Eq. (3) and achieved successful predictions.

Adjustment of the relative position between receptor and ligand

The interface of the docked scaffold protein was designed to bind the receptor protein, and binding affinity (ΔG1) was calculated. The ligand protein is subsequently subjected to a series of rotations and translations. These displacements are performed around a randomly selected axis through the geometry center of the ligand with a step size of (−5° < а <5°) and d (−1.0 Å < d < 1.0 Å). The interface of the ligand protein was designed after the adjustment. A binding affinity (ΔG2) is then calculated. If ΔΔG = ΔG2 − ΔG1 < 0, the move was accepted. Otherwise, the move was accepted with a probability estimated by the following expression: exp[–ΔΔG/T]. The annealing temperature T was set to 2 initially and scaled by 0.98 after each move. The procedure was repeated 100 times. It should be mentioned that the mutations at each move were not transferred to the next step. The non-interface residues of the scaffold were always the same as the native protein.

Flexible backbone design

We selected 50 conformers from the MD simulation of the complex. For each conformer, we repeated to design the interface five times. The one calculated as the lowest binding affinity was selected from the 50 × 5 designed complexes. We also performed side-chain modeling of interface residues five times for each conformer. The side-chain conformation was modeled with rotamers as in the design procedure but not mutated. The one calculated as the lowest binding affinity was compared with that from the interface design.

Protein structures

As in the previous study (Liang et al., 2009), three protein–protein complexes were considered in this study: barnase–barstar (PDB code: 1brs) from the RNAase-inhibitor family; hen egg lysozyme–Fv fragment of antibody D1.3 (PDB code: 1vfb) from antigen-antibody family; and trypsin–bovine pancreatic trypsin inhibitor (PDB code: 2ptc) from the proteinase-inhibitor family. The interface of barstar, antibody D1.3 and BPTI was designed in the conformers from the MD simulation of the native structure. We also designed binding site in the scaffold proteins to bind the unbound structure of barnase (1a2p), trypsin (2ptn) or lysozyme (8lyz). Three scaffold proteins tenascin (1ten), CheY (1tmy) and DNA-binding domain of MBP1 (1bm8) were used in this study. 1ten is a β-protein, whereas most secondary structure components are helixes on the surface of 1tmy and 1bm8 with only 2–3 solvent-accessible β-strands.

MD simulations

The following protocol for setting up and running the MD simulations was followed in all cases. Hydrogen atoms were added to the protein using the PROTONATE program, which is part of the AMBER 9 suite of programs. The AMBER force field parameters were assigned to all atoms using the ‘parm99’ set of parameters. The SYBYL 8.0 program was used for the manipulation and visualization of all structures and for protonation of the bound ligands. The program LEaP was used to neutralize the complexes. The complexes of barstar–barnase, antibody D1.3–lysozyme and BPTI–trypsin were immersed in a box of TIP3P water molecules such that no atom in the complex was within 12 Å for all complexes from any side of the box. All bonds involving hydrogen atoms were constrained by using the SHAKE algorithm, affording the use of a 2 fs time-step. The particle mesh Ewald method was used to treat long-range electrostatics. Water molecules were first energy minimized and equilibrated by running a short simulation with the complex fixed by using Cartesian restraints. This was followed by a series of energy minimizations in which the Cartesian restraints were gradually relaxed from 500 kcal/Å, and the system was subsequently gradually heated to 300 K via a 48 ps MD run. Another 1 ns simulation was carried out at 300 K for further equilibration.

MD-based free energy calculations

The method for determining the binding free energy following the MM-PBSA approach has been described in the past (Srinivasan et al., 1998; Massova and Kollman, 2000). The binding free energy is expressed as:  

(4)
formula
where 
(5)
formula
where 〈EINT〉 is the internal potential energy determined using the AMBER force field, 〈Gsolvnp〉 the non-polar contribution to the solvation free energy and 〈Gsolvele〉 the electrostatic contribution to the free energy and 〈Stot〉 the configurational entropy of the system given by  
(6)
formula
where 〈Strans〉, 〈Srot〉 and 〈Svib〉 are the translational, rotational and vibrational components of the absolute entropy. In this work, the MM-PBSA free energies are determined by extracting 320 snapshots from the multiple trajectories for each complex, the same number of snapshots was used to compute the entropy change with normal mode analysis.

The electrostatic contribution to the solvation free energy is determined by using the PBSA program in AMBER 9 which numerically solves the Poisson–Boltzmann equations to determine the electrostatic contribution to the solvation free energy. A 0.5 Å grid size was used, and the dielectric constant for the solute and solvent was set to 1 and 80, respectively. The optimized atomic radii set in AMBER 9 were used, and partial charges were taken from Cornell et al. for standard residues. The non-polar contribution to the solvation free energy was determined by using the MOLSURF program, which is also part of the AMBER 9 suite of programs. The entropy was computed for each snapshot from normal mode analyses by using the NMODE program within the AMBER package.

Results and discussion

Performance of scoring function in CAPRI

CAPRI is a community wide experiment to assess protein docking methods and protein–protein interaction energy functions (http://capri.ebi.ac.uk/). Groups participate in CAPRI in two capacities: predictor or scorer. As a predictor users submit blind structure predictions for protein–protein target complexes based on unbound structures; as a scorer CAPRI users select near-native decoys from a pool of decoys submitted by predictors. We participated in CAPRI rounds 10, 11 and 12 both as predictor and scorer using Eq. (3). The following steps were used to generate 10 top models for the prediction and scoring tests.

In the prediction test, we first collect all possible biological and structural information from the literature concerning the possible binding site. Then ZDOCK is used to generate top 2000 decoys, and the hypothesized non-binding regions are blocked (ZDOCK applies special penalty for the non-binding or blocked regions). Default parameters for the program were used for the docking, except for the degree of rotational sampling intervals, which was set to 10. The interface residues of the 2000 complexes are optimized by side-chain modeling, and the complexes are further optimized by 50 steps of ABNR energy minimization (see ‘Methods’ section). Finally, the binding affinity was evaluated by Eq. (3).

We also participated in CAPRI as a scorer. Here, the decoys to test scoring functions are provided by CAPRI participants and refined by various methods. We found that our scoring function is sensitive to the structural refinement methods (Liang et al., 2007). Energy minimization resulted in lower free energy. The computed binding affinities are not comparable for these decoys. We thus superimpose the unbound monomers to the corresponding protein in the complex. The complexes of the superimposed unbound monomers are then refined and evaluated as we did in the prediction test.

The incorporation of existing experimental data during computational prediction may significantly improve predictions in CAPRI competitions. For example, in the case of target 26 TolB/Pal in round 10, it is well known from X-ray diffraction and site-directed mutagenesis data (Ray et al., 2000) that the C-terminal β-propeller structure of TolB is responsible for the interaction with Pal. The region of Pal that interacts with TolB is delimited by residues 89–104 and 126–130 (Clavel et al., 1998). Hence, regions that are not involved in the TolB/Pal complex formation are blocked during the ZDOCK sampling. In the case of target 28 E2-25K/Ubc9 in round 11 (target 27 could be contact of E2-25K/Ubc9 in the crystal structure and none of the participants obtained a successful prediction), Ubc9 is required for SUMO (small ubiquitin-related modifier) modification of the ubiquitin-conjugating enzyme E2-25K (Johnson, 2004). The acceptor lysine for SUMO conjugation in E2-25K is Lys14 (Pichler et al., 2005). The donor ubiquitin is attached to Cys93 in Ubc9. We guess that Lys14 of E2-25K should be close to Cys93 of Ubc9 in the complex of E2-25K/Ubc9. Thus, residues 1–59 and 99–123 in Ubc9, which are on the opposite side of Cys93 in the crystal structure, are blocked before docking sampling. Similarly, residues 78–98 and 121–199 in E2-25K are blocked. One unbound structure of the target 29 in round 12 was derived from homology modeling. None of the participants obtained a successful prediction, and the results are not discussed. The overall performance of our best predictions for targets 26 and 28 is summarized in Table I. There are typically 40 groups that participate in CAPRI in each round. We were ranked number 7 and number 3 for target 26 and target 28, respectively, in term of global root mean squared (r.m.s.) deviation. For both targets, the docking accuracy of the prediction test was found to be similar to that of the scoring test. It is encouraging that the scoring accuracy was similar to that of docking accuracy given that no experimental data were used in the scoring.

Table I

The best model according to the overall r.m.s. deviation value of CAPRI target 26 and target 28

Target Test F (nat)a f (non-nat)b L (r.m.s.d.)c I (r.m.s.d.)d Classification Ranke 
T26 Prediction 0.489 0.281 5.430 1.862 Medium 
Scoring 0.511 0.294 4.614 2.225 Medium 
T28 Prediction 0.463 0.694 5.644 2.927 Acceptable 
Scoring 0.366 0.681 10.405 2.804 Acceptable 11 
Target Test F (nat)a f (non-nat)b L (r.m.s.d.)c I (r.m.s.d.)d Classification Ranke 
T26 Prediction 0.489 0.281 5.430 1.862 Medium 
Scoring 0.511 0.294 4.614 2.225 Medium 
T28 Prediction 0.463 0.694 5.644 2.927 Acceptable 
Scoring 0.366 0.681 10.405 2.804 Acceptable 11 

aThe fraction of ligand–receptor contacts that is also found in the native structure.

bThe fraction of ligand–receptor contacts that is found, but that is not present in the native structure.

cThe r.m.s. deviation between target and predicted ligand molecule(s) after the receptor molecules are fitted.

dThe r.m.s. deviation value of the interface backbone atoms.

eThe rank in the CAPRI evaluation team.

Alanine scanning at the binding interface

In contrast to computing absolute binding affinities, it is less challenging to determine the binding free energy change due to mutation. The interaction energy between two proteins was calculated similar to that between the position rotamer and other parts of protein in side-chain modeling. We assume that the mutations do not cause any structural alteration of the protein complex. The binding energy of the complex mutant is determined using the crystal structure of wild type except that the predicted residue is replaced by alanine.

We assessed the performance of our scoring function in predicting the binding free energy changes upon alanine mutagenesis for barnase–barstar and lysozyme–antibody D1.3 complexes (Figs 1 and 2). Experimental values were obtained from previous studies of these extensively studied systems (Schreiber and Fersht, 1995; Dall'Acqua et al., 1998). A linear fit yielded a correlation coefficient of 0.84 between the computed and experimental binding free energy changes in the case of lysozyme–antibody D1.3. We obtained a lower correlation coefficient (0.72) for barnase–barstar. This slightly lower performance may be due to the fact that the electrostatic interactions for core and surface residues were treated equally in our scoring function. The interaction energy of a buried salt bridge in the barnase–barstar complex may therefore have been underestimated. If the four mutations (charged residue→Ala) at the relatively buried position are excluded, we obtain a correlation efficient of 0.97 for the other six mutations. The prediction accuracy is similar to that of Kortemme and Baker (2002). They obtained a correlation coefficient of 0.83 for lysozyme–antibody D1.3 and 0.96 for the interface residues of barnase–barstar (Kortemme and Baker, 2002). It is of interest to note that these investigators obtained good correlation when solvent-exposed salt bridges are excluded, which is in contrast to our results, where better correlations were obtained when buried salt bridges were not considered.

Fig. 1

Comparison of calculated interaction energy score and experimental binding free-energy changes upon alanine mutation. (A) Lysozyme–antibody D1.3 and (B) barnase–barstar. The regression lines were (A) ΔEcalc = 2.34 ΔGexp − 0.012 and (B) ΔEcalc = 1.04 ΔGexp + 1.05.

Fig. 1

Comparison of calculated interaction energy score and experimental binding free-energy changes upon alanine mutation. (A) Lysozyme–antibody D1.3 and (B) barnase–barstar. The regression lines were (A) ΔEcalc = 2.34 ΔGexp − 0.012 and (B) ΔEcalc = 1.04 ΔGexp + 1.05.

Fig. 2

Comparison of native and redesigned lysozyme–antibody D1.3 interface. Proteins are shown in deep teal and orange ribbon representations for the receptor and ligand. The upper and lower panels show the native and mutated interfaces, respectively. Several residues at the interface are shown in capped stick representation color-coded according to atom types (yellow, red and blue for C, O and N, respectively). A colour version of this figure is available as supplementary data at PEDS online.

Fig. 2

Comparison of native and redesigned lysozyme–antibody D1.3 interface. Proteins are shown in deep teal and orange ribbon representations for the receptor and ligand. The upper and lower panels show the native and mutated interfaces, respectively. Several residues at the interface are shown in capped stick representation color-coded according to atom types (yellow, red and blue for C, O and N, respectively). A colour version of this figure is available as supplementary data at PEDS online.

Interface redesign with flexible backbone

Before we embark on the grafting of protein interfaces, we sought to put our algorithm and scoring function to the test by redesigning the interfaces of three protein complexes, barnase–barstar (PDB ID: 1BRS), trypsin–BPTI (PDB ID: 2PTC) and lysozyme–antibody D1.3 (PDB ID: 1VFB). The process consists of systematically mutating a number of residues at the interface of these complexes until a more stable complex is created. The algorithm is described in detail in the ‘Methods’ section.

In a previous study, we had redesigned the interface of a number of protein–protein complexes (Liang et al., 2009). The mutations were performed using a Monte Carlo-based algorithm with the binding affinity assessed by our scoring function described in the ‘Methods’ section. A number of complexes that exhibited greater affinity than the native structure were identified. Interestingly, the designer interfaces showed high sequence similarity with the native complex. Despite the successes, we showed in the same study that the lowest energy complexes that emerged from our design algorithm did not always result in stable complexes when subjected to explicit-solvent MD simulations. In fact, a number of these complexes completely dissociated during the trajectory (Liang et al., 2009).

The design process from our previous work was performed on a fixed receptor structure. Here, we extend on this work by explicitly incorporating backbone flexibility of the receptor with our design protocol. We start with a simple redesign of the interface of the three complexes, namely barnase–barstar (PDB ID: 1BRS), trypsin–BPTI (PDB ID: 2PTC) and lysozyme–antibody D1.3 (PDB ID: 1VFB). The design algorithm follows closely the one that we have employed previously, except that it is performed on 50 conformers that are collected from a 6 × 5 ns explicit-solvent MD simulation of the complexes. The designed complexes with the highest binding affinity were compared with the conformer with the highest binding affinity after side-chain modeling of the interface residues.

An example of a stable interface that emerged from redesigning the barnase–barstar, trypsin–BPTI and lysozyme–antibody D1.3 is shown in Fig. 2. For the barnase (1brs) and trypsin (2ptc) complexes, the designed sequences showed greater affinity over the native complex by nearly 2 kcal/mol as determined by Eq. (3) (Table II). The small improvement in affinity is not surprising given the already tight interactions between the binding partners of these complexes. In fact, these interactions are among the tightest that are observed in nature. On the other hand, we obtained a greater improvement for those complexes that were designed starting with the lysozyme–antibody D1.3 complex as shown in Fig. 2A; the most stable complex showed a 4.4 kcal/mol increase in affinity over the native structure (Fig. 2B).

Table II

Comparison of binding free energy between native and designed sequences in conformers from the MD simulation of the native complex

Complex PDB ID Modeling Design Difference 
Barnase–barstar 1BRS −16.1 −18.0 −1.9 
Trypsin–BPTI 2PTC −18.8 −19.2 −0.4 
Lysozyme–antibody D1.3 1VFB −9.9 −14.3 −4.4 
Complex PDB ID Modeling Design Difference 
Barnase–barstar 1BRS −16.1 −18.0 −1.9 
Trypsin–BPTI 2PTC −18.8 −19.2 −0.4 
Lysozyme–antibody D1.3 1VFB −9.9 −14.3 −4.4 

Values are reported in kcal/mol.

To probe for their stability in solution, designer complexes are subjected to extensive explicit-solvent MD simulations for 30 ns each (six trajectories, 5 ns in length). The motion experienced by these complexes was monitored by following the evolution of the r.m.s. deviation of the backbone atoms over the course of the trajectories as shown in Fig. 3. It was highly encouraging that the r.m.s. deviation of the overall complex in all three complexes with redesigned interfaces remained within 2 Å, suggesting that they are likely to be stable in solution.

Fig. 3

Structural changes of the barnase-barstar (solid), lysozyme-antibody D1.3 (dash), and trypsin-BPTI (dotted) complexes with redesigned interfaces. Root mean squared deviation is monitored over the course of six 5 ns molecular dynamics simulations.

Fig. 3

Structural changes of the barnase-barstar (solid), lysozyme-antibody D1.3 (dash), and trypsin-BPTI (dotted) complexes with redesigned interfaces. Root mean squared deviation is monitored over the course of six 5 ns molecular dynamics simulations.

Regardless of the structural integrity that was exhibited by the designer complexes during the MD simulations, we sought to gain a deeper understanding of the forces that led to stable complexes by performing extensive free energy calculations and comparing the results to the native structure. These free energy calculations follow the MM-PBSA approach whereby the electrostatic, non-polar and entropy components are averaged over multiple conformers collected from the MD simulation (Srinivasan et al., 1998; Massova and Kollman, 2000). The results for the three complexes are shown in Table III. In the case of the system with redesigned barnase–barstar interface, the overall free energy is significantly greater in the designer complex than the native complex (10.7 versus 2.9 kcal/mol). This is despite the fact that both native and designer complexes are stable during the MD simulation. A close inspection of the components of the free energy shows that the electrostatic ΔGPBELE is primarily responsible for the better overall binding free energy of the native complex, whereas entropy opposes binding. Interestingly, the non-polar component of the free energy ΔGNP is similar in both cases, as evidenced by comparable contributions from non-polar solvation (ΔGSA) and potential energy from non-bonded interactions (ΔGVDW). This observation is consistent with our previous finding that non-polar interactions were a strong predictor of stability in solution, even more than the overall free energy or its other components (Liang et al., 2009).

Table III

Free energy calculations for designer complexes and complexes with transferred interfacea

Receptor PDB Type Stability ΔGELE ΔGPB ΔGPBELE ΔGSA ΔGVDW ΔGNP TΔSNM ΔGcalc 
Barnase – Grafted STABLE −312.2 ± 1.5 385.9 ± 1.4 73.7 ± 0.1 −10.7 ± 0.0 −94.9 ± 0.3 −105.6 ± 0.3 −42.6 ± 0.4 10.7 ± 0.9 
Barnase 1BRS Native Stable −527.0 ± 2.5 588.6 ± 2.3 61.6 ± 0.9 −12.2 ± 0.03 −93.5 ± 0.3 −105.7 ± 0.3 −45.6 ± 1.2 2.9 ± 1.4 
Lysozyme – Grafted Stable 299.1 ± 2.2 −226.4 ± 2.3 72.7 ± 0.1 −12.3 ± 0.1 −77.0 ± 0.6 −89.3 ± 0.6 −46.0 ± 0.6 29.4 ± 1.1 
Lysozyme 1VFB Native Stable 121.1 ± 1.1 −68.6 ± 1.0 52.5 ± 0.5 −10.8 ± 0.04 −73.7 ± 0.3 −84.5 ± 0.3 −43.5 ± 1.3 11.5 ± 1.4 
Trypsin – Grafted Stable 183.9 ± 1.1 −177.9 ± 1.0 6.1 ± 0.1 −9.3 ± 0.0 −95.7 ± 0.4 −105.0 ± 0.4 −35.8 ± 0.4 −63.1 ± 0.7 
Trypsin 2PTC Native Stable 234.3 ± 1.5 −224.0 ± 1.3 10.3 ± 0.6 −11.1 ± 0.03 −90.4 ± 0.3 −101.5 ± 0.3 −47.5 ± 1.3 −43.7 ± 1.4 
Receptor PDB Type Stability ΔGELE ΔGPB ΔGPBELE ΔGSA ΔGVDW ΔGNP TΔSNM ΔGcalc 
Barnase – Grafted STABLE −312.2 ± 1.5 385.9 ± 1.4 73.7 ± 0.1 −10.7 ± 0.0 −94.9 ± 0.3 −105.6 ± 0.3 −42.6 ± 0.4 10.7 ± 0.9 
Barnase 1BRS Native Stable −527.0 ± 2.5 588.6 ± 2.3 61.6 ± 0.9 −12.2 ± 0.03 −93.5 ± 0.3 −105.7 ± 0.3 −45.6 ± 1.2 2.9 ± 1.4 
Lysozyme – Grafted Stable 299.1 ± 2.2 −226.4 ± 2.3 72.7 ± 0.1 −12.3 ± 0.1 −77.0 ± 0.6 −89.3 ± 0.6 −46.0 ± 0.6 29.4 ± 1.1 
Lysozyme 1VFB Native Stable 121.1 ± 1.1 −68.6 ± 1.0 52.5 ± 0.5 −10.8 ± 0.04 −73.7 ± 0.3 −84.5 ± 0.3 −43.5 ± 1.3 11.5 ± 1.4 
Trypsin – Grafted Stable 183.9 ± 1.1 −177.9 ± 1.0 6.1 ± 0.1 −9.3 ± 0.0 −95.7 ± 0.4 −105.0 ± 0.4 −35.8 ± 0.4 −63.1 ± 0.7 
Trypsin 2PTC Native Stable 234.3 ± 1.5 −224.0 ± 1.3 10.3 ± 0.6 −11.1 ± 0.03 −90.4 ± 0.3 −101.5 ± 0.3 −47.5 ± 1.3 −43.7 ± 1.4 

aUnits in kcal/mol; ΔGELE, non-solvent electrostatic potential energy; ΔGPB, electrostatic contributions to the solvation free energy calculated with the Poisson–Boltzmann equation; ΔGPBELE, electrostatic potential energy; ΔGSA, non-polar contributions to solvation free energy; ΔGVDW, van der Waals potential energy; ΔGNP, non-polar potential energy; TΔSNM, the entropic contribution calculated with normal mode analysis to the free energy of binding; ΔGcalc, calculated binding free energy.

In the case of the lysozyme–antibody D1.3 designer complexes, the free energy of the grafted complex (ΔGcalc = 29.4 kcal/mol) was also less favorable than the native complex (ΔGcalc = 11.5 kcal/mol). Electrostatics and entropy are mainly responsible for this, as these components are strongly more favorable in the native complex. In contrast, the non-polar component (ΔGNP) was more favorable in the designer complex by nearly 5 kcal/mol. Unlike the case of barnase–barstar, the more favorable non-polar component was due to both the non-bonded potential interaction energy (ΔGVDW) and the larger burial of solvent-accessible surface area (ΔGSA) upon complex formation. This again shows that the non-polar term of the free energy (ΔGNP) is a strong predictor of stability.

Finally, the free energy data for the flexible backbone redesign of the trypsin–BPTI complex interface resulted in a complex with a binding affinity that was significantly more favorable than the native complex. This is in contrast with the barnase and lysozyme designer complexes that exhibited lower affinity than the native complex. Interestingly, the data for the trypsin designer complex show that every component of the free energy is more favorable. The electrostatic energy ΔGPBELE was nearly 4 kcal/mol lower than the native complex. This primarily was due to the significantly more favorable Coulomb energy (ΔGELE), which compensated for the unfavorable contributions from the electrostatic component of the solvation free energy ΔGPB. The non-polar component, which includes the van der Waals energy as well as the non-polar solvation free energy, was more favorable for the designer complex, an observation that has been made so far for every complex that remained stable during MD. Finally, the complex even exhibited more favorable entropy of binding, suggesting that less tightening of the structure occurs during binding.

Interface design of unrelated receptors

Although the redesign of existing interfaces is valuable and has been successfully accomplished in the past by a number of investigators (Ogata et al., 2003; Ali et al., 2005; Green et al., 2006; Lippow et al., 2007; Yosef et al., 2009), the redesign of the surface of a protein to bind to evolutionarily unrelated proteins is a significantly greater challenge with only a handful of successes reported in the literature (Liu et al., 2007). In a previous study, we proposed a design algorithm that was primarily guided by transferring a binding epitope in the ligand of the native complex to a stably folded scaffold protein. The process consisted of first defining a set of hot-spot residues or binding epitope on the ligand and then scanning the surface of a large set of scaffolds for a set of residues that adopt the arrangement of the hot-spot residues in the native complex. We identified a number of candidates and re-grafted their surface such that they would bind to barnase, trypsin or lysozyme (Liang et al., 2009). The resulting complexes were all subjected to extensive explicit-solvent MD simulations, and we found a number of complexes were stable, whereas others dissociated during the trajectories.

A limitation of the aforementioned approach is that the X-ray or NMR structure of the receptor in complex with a ligand is prerequisite. In many situations, such a structure of the complex is not available and only that of the receptor can be accessed. In this study, we address this limitation by performing our computational design using strictly the structure of the receptor. This was primarily accomplished by unrestrained molecular docking of the scaffold to the receptor, with the single exception that the docking region on the receptor is confined to a smaller area of the surface that typically corresponds to the known binding site. Second, we follow a flexible backbone approach whereby multiple conformers of the receptor are employed during the design.

In our previous work (Liang et al., 2009), we had found that a number of designer complexes were stable during MD; these included the complex between lysozyme and grafted MPB1 (1VFB_1BM8), trypsin and grafted CheY (2PTC_1TMY) or tenascin (2PTC_1TEN). None of the designer complexes for barnase (1BRS) resulted in stable complexes. Here, we employ our flexible backbone approach and attempt to re-graft MPB1, CheY and tenascin to bind barnase, lysozyme and trypsin resulting in 3 × 3 = 9 designer complexes. The process involved first the docking of the ligands to their receptors, followed by a Monte Carlo-based algorithm to redesign their interface through a series of mutations carried out on multiple residues simultaneously.

Typically, it was possible to identify—among the top 10 docked candidates—at least one complex that did not exhibit significant gaps at the interface. However, a number of side chain–backbone and side chain–side chain clashes were found. This is not unexpected since ZDOCK2.3 was designed for docking of unbound structures that follows a rigid docking approach. These clashes were eventually eliminated during the design process. It was observed that the binding free energy of the docking was positive using our all-atom scoring function, but was found to decrease after mutation of interface residues and refinement of the relative position between receptor and ligand. Clearly, the rigid-body transformation improves the binding affinity of the designed scaffold proteins as shown in Table IV. The calculated binding affinity was comparable with that of native complexes. Interestingly, Glu35 on the designed scaffold of 1tmy that was grafted to bind barnase is located in a similar position as Asp39 of barstar, which is a well-known hot-spot residue in the barnase–barstar complex with a binding free energy change >5 kcal/mol upon mutation to Ala. Lys70 of the designed scaffold of MBP1 (1bm8) to bind trypsin is located at the same position as Lys15, a hot-spot residue for the trypsin–BPTI complex.

Table IV

Comparison of binding affinity for the designed complexes before and after adjustment of relative position

Scaffold Ligand
 
 Tenascin (1ten)
 
CheY (1tmy)
 
MBP1 (1bm8)
 
 Beforea Afterb Before After Before After 
Barnase (1a2p) −0.3 −7.9 −9.9 −13.8 −2.6 −14.5 
Trypsin (2ptn) −6.3 −14.3 −6.6 −13.4 3.0 −9.3 
Lysozyme (8lyz) −2.3 −9.8 −1.9 −9.0 0.3 −9.4 
Scaffold Ligand
 
 Tenascin (1ten)
 
CheY (1tmy)
 
MBP1 (1bm8)
 
 Beforea Afterb Before After Before After 
Barnase (1a2p) −0.3 −7.9 −9.9 −13.8 −2.6 −14.5 
Trypsin (2ptn) −6.3 −14.3 −6.6 −13.4 3.0 −9.3 
Lysozyme (8lyz) −2.3 −9.8 −1.9 −9.0 0.3 −9.4 

aThe binding free energy was calculated when the interface was designed in the docking complex.

bThe binding free energy was calculated after the relative position between receptor and ligand proteins were adjusted by Monte Carlo methods.

To identify those designer complexes that are likely going to be stable in solution, we subjected all nine designer complexes (shown in Fig. 4) to explicit-solvent MD simulations for 30 ns each, for a total of nearly 270 ns. The structural integrity of these complexes is monitored over the course of the trajectories (Fig. 5) and conformers sampled during the trajectories are shown in Fig. 6. Interestingly, some of the designer complexes were stable, whereas others dissociated during the MD simulations as evidenced by the significant increase in the r.m.s. deviation away from the initial complex as monitored with respect to time (Fig. 5). In the case of barnase, only one of the ligands showed stability toward it as shown in Fig. 6A (top). This is highly interesting since in our last effort (Liang et al., 2009), none of the designer complexes with barnase were stable. This shows that the docking approach may in some cases be more appropriate. An example of an unstable complex involving barnase is shown in Fig. 6A (bottom). The complexes involving lysozyme resulted in two stable complexes and one unstable complex, as illustrated in Fig. 6B (top and bottom, respectively). Grafted MPB1 resulted in a stable complex with lysozyme, consistent with the results from our previous study following a hot-spot driven design algorithm. Finally, two of the complexes with trypsin dissociated (tenascin and CheY). It is interesting to note that in our previous effort, designer complexes of trypsin with both CheY (1tmy) and tenascin (1ten) resulted in stable complexes. It is possible that in this particular case, our previous hot-spot driven algorithm was more appropriate. The docking approach may not have efficiently exploited those residues that are critical for binding in trypsin.

Fig. 4

Stereoview of the three dimensional structure of designer complexes that involve barnase, lysozyme and trypsin as receptors and MBP1, tenascin and CheY as ligands. (A) 1BRS_1BM8, (B) 1BRS_1TEN, (C) 1BRS_1TMY, (D) 1VFT_1BM8, (E) 1VFT_1TEN, (F) 1VFB_1TMY, (G) 2PTC_1BM8, (H) 2PTC_1TEN and (I) 2PTC_1TMY. In each case, receptor is shown in deep teal and ligand in orange ribbon rendering. A colour version of this figure is available as supplementary data at PEDS online.

Fig. 4

Stereoview of the three dimensional structure of designer complexes that involve barnase, lysozyme and trypsin as receptors and MBP1, tenascin and CheY as ligands. (A) 1BRS_1BM8, (B) 1BRS_1TEN, (C) 1BRS_1TMY, (D) 1VFT_1BM8, (E) 1VFT_1TEN, (F) 1VFB_1TMY, (G) 2PTC_1BM8, (H) 2PTC_1TEN and (I) 2PTC_1TMY. In each case, receptor is shown in deep teal and ligand in orange ribbon rendering. A colour version of this figure is available as supplementary data at PEDS online.

Fig. 5

Mean r.m.s. deviation over six trajectories for (A) barnase, (B) lysozyme, and (C) trypsin designer complexes. The lines are shown in solid, dash and dotted corresponding to MBP1, tenascin, and CheY.

Fig. 5

Mean r.m.s. deviation over six trajectories for (A) barnase, (B) lysozyme, and (C) trypsin designer complexes. The lines are shown in solid, dash and dotted corresponding to MBP1, tenascin, and CheY.

Fig. 6

Stereoview of the three-dimensional structures from MD simulations of stable (top) and unstable (bottom) complexes for (A) barnase, (B) lysozyme and (C) trypsin. A colour version of this figure is available as supplementary data at PEDS online.

Fig. 6

Stereoview of the three-dimensional structures from MD simulations of stable (top) and unstable (bottom) complexes for (A) barnase, (B) lysozyme and (C) trypsin. A colour version of this figure is available as supplementary data at PEDS online.

To gain a better understanding of the potential reasons for instability in some of these complexes, we collected structures of the nine complexes after 1 ns of equilibration and performed a detailed analysis of the binding interface using the protein–protein interaction server (http://www.bioinformatics.sussex.ac.uk/protorp/). The data are provided in Table V. It is generally assumed that greater gaps at the interface should lead to less stable complexes, but this was not supported by our data. For example, the trypsin–MBP1 (2PTC_1BM8) complex was found to be stable, but it contained more gaps than any of the other systems. Conversely, the complex between barnase and MBP1 (1BRS_1BM8) exhibited the smallest number of gaps, but was unstable. When the size of the binding interface is taken into consideration by taking the ratio of gap volume over the interaction surface, there is no trend that emerges. Gaps at protein–protein binding interface provide space for water molecules to occupy and potentially act as a bridge between the binding partners. Here, more water molecules at the interface did not necessarily result in stable complexes. For example, the trypsin–MBP1 (2PTC_1BM8) complex showed only 5 bridging waters and was stable, whereas the trypsin–tenascin (2PTC_1TEN) complex had 19 bridging waters but was unstable. Finally, the number of salt bridges or hydrogen bonds was not a good predictor of stability, as the unstable lysozyme–tenascin (1VFB_1TEN) complex showed 30 salt bridges and 12 hydrogen bonds, but the stable trypsin–MBP1 (2PTC_1BM8) showed only 11 salt bridges and 9 hydrogen bonds.

Table V

Protein–protein interface properties for designer complexes

Receptor Ligand Stability Vgap3)a VIgap (Å)a Bridging water Salt bridge H-bond 
Barnase MBP1 Unstable 3152.3 1.56 10 
Barnase Tenascin Stable 4029.8 2.18 14 24 23 
Barnase CheY Unstable 3948.8 2.51 15 
Lysozyme MBP1 Unstable 3199.5 2.50 11 
Lysozyme Tenascin Unstable 3547.1 2.67 30 12 
Lysozyme CheY Stable 2669.6 1.99 10 
Trypsin MBP1 Stable 5248.1 3.62 11 
Trypsin Tenascin Unstable 3543.8 1.80 19 
Trypsin CheY Unstable 4077.0 2.87 16 
Receptor Ligand Stability Vgap3)a VIgap (Å)a Bridging water Salt bridge H-bond 
Barnase MBP1 Unstable 3152.3 1.56 10 
Barnase Tenascin Stable 4029.8 2.18 14 24 23 
Barnase CheY Unstable 3948.8 2.51 15 
Lysozyme MBP1 Unstable 3199.5 2.50 11 
Lysozyme Tenascin Unstable 3547.1 2.67 30 12 
Lysozyme CheY Stable 2669.6 1.99 10 
Trypsin MBP1 Stable 5248.1 3.62 11 
Trypsin Tenascin Unstable 3543.8 1.80 19 
Trypsin CheY Unstable 4077.0 2.87 16 

aVgap, gap volume between two proteins; VIgap, gap volume index, defined as the gap volume over the interface solvent-accessible surface area.

To gain a better understanding of the forces that lead to stable complexes, the components of the free energy for each of the nine complexes were computed as shown in Table VI. The designer complexes for the barnase receptor resulted in two unstable and one stable complex. The stable complex (1BRS_1TEN) showed an overall binding affinity (ΔGcalc) that is significantly more favorable than the unstable complexes (−11 versus 34 and 41, respectively). The non-polar component, which is usually predictive of stability, was, not surprisingly, more favorable in the stable complex. The electrostatic component of the free energy ΔGPBELE was also significantly more favorable in the stable complex. However, it appears that the stable complex incurred greater entropic penalty for binding.

Table VI

Free energy calculations for native and designer complexesa

Receptor PDB Stability ΔGELE ΔGPB ΔGPBELE ΔGSA ΔGVDW ΔGNP T ΔSNM ΔGcalc 
Barnase 1BRS_1BM8 Unstable 94.7 ± 2.2 −2.7 ± 2.5 92.1 ± 0.1 −11.3 ± 0.1 −92.8 ± 0.9 −104.1 ± 0.9 −46.2 ± 0.4 34.2 ± 1.5 
Barnase 1BRS_1TEN Stable −577.3 ± 2.3 621.5 ± 2.4 44.2 ± 0.1 −12.2 ± 0.0 −93.0 ± 0.4 −105.2 ± 0.4 −49.6 ± 0.4 −11.4 ± 0.9 
Barnase 1BRS_1TMY Unstable −255.8 ± 1.9 338.0 ± 2.3 82.2 ± 0.1 −10.5 ± 0.0 −74.0 ± 0.5 −84.6 ± 0.5 −43.6 ± 0.3 41.2 ± 0.8 
Lysozyme 1VFB_1BM8 Unstable 628.5 ± 2.5 −561.9 ± 2.4 66.6 ± 0.1 −8.2 ± 0.0 −61.8 ± 0.4 −69.9 ± 0.4 −33.7 ± 0.4 30.3 ± 0.7 
Lysozyme 1VFB_1TEN Unstable −1240.1 ± 2.3 1271.8 ± 2.3 31.6 ± 0.1 −7.2 ± 0.0 −50.0 ± 0.4 −57.1 ± 0.4 −33.6 ± 0.4 8.1 ± 0.8 
Lysozyme 1VFB_1TMY Stable −133.9 ± 2.0 179.8 ± 2.1 45.9 ± 0.1 −8.2 ± 0.0 −64.8 ± 0.3 −73.0 ± 0.3 −32.8 ± 0.3 5.7 ± 0.7 
Trypsin 2PTC_1BM8 Stable 267.4 ± 2.2 −222.3 ± 2.1 45.1 ± 0.1 −10.1 ± 0.0 −92.0 ± 0.3 −102.1 ± 0.3 −40.2 ± 0.5 −16.8 ± 0.8 
Trypsin 2PTC_1TEN Unstable −597.1 ± 5.3 664.8 ± 7.0 67.7 ± 0.4 −9.2 ± 0.2 −66.2 ± 1.8 −75.4 ± 1.8 −43.5 ± 0.9 35.8 ± 1.2 
Trypsin 2PTC_1TMY Unstable −157.7 ± 3.3 230.9 ± 3.6 73.2 ± 0.2 −8.6 ± 0.1 −61.6 ± 0.5 −70.1 ± 0.5 −39.7 ± 0.6 42.8 ± 1.1 
Receptor PDB Stability ΔGELE ΔGPB ΔGPBELE ΔGSA ΔGVDW ΔGNP T ΔSNM ΔGcalc 
Barnase 1BRS_1BM8 Unstable 94.7 ± 2.2 −2.7 ± 2.5 92.1 ± 0.1 −11.3 ± 0.1 −92.8 ± 0.9 −104.1 ± 0.9 −46.2 ± 0.4 34.2 ± 1.5 
Barnase 1BRS_1TEN Stable −577.3 ± 2.3 621.5 ± 2.4 44.2 ± 0.1 −12.2 ± 0.0 −93.0 ± 0.4 −105.2 ± 0.4 −49.6 ± 0.4 −11.4 ± 0.9 
Barnase 1BRS_1TMY Unstable −255.8 ± 1.9 338.0 ± 2.3 82.2 ± 0.1 −10.5 ± 0.0 −74.0 ± 0.5 −84.6 ± 0.5 −43.6 ± 0.3 41.2 ± 0.8 
Lysozyme 1VFB_1BM8 Unstable 628.5 ± 2.5 −561.9 ± 2.4 66.6 ± 0.1 −8.2 ± 0.0 −61.8 ± 0.4 −69.9 ± 0.4 −33.7 ± 0.4 30.3 ± 0.7 
Lysozyme 1VFB_1TEN Unstable −1240.1 ± 2.3 1271.8 ± 2.3 31.6 ± 0.1 −7.2 ± 0.0 −50.0 ± 0.4 −57.1 ± 0.4 −33.6 ± 0.4 8.1 ± 0.8 
Lysozyme 1VFB_1TMY Stable −133.9 ± 2.0 179.8 ± 2.1 45.9 ± 0.1 −8.2 ± 0.0 −64.8 ± 0.3 −73.0 ± 0.3 −32.8 ± 0.3 5.7 ± 0.7 
Trypsin 2PTC_1BM8 Stable 267.4 ± 2.2 −222.3 ± 2.1 45.1 ± 0.1 −10.1 ± 0.0 −92.0 ± 0.3 −102.1 ± 0.3 −40.2 ± 0.5 −16.8 ± 0.8 
Trypsin 2PTC_1TEN Unstable −597.1 ± 5.3 664.8 ± 7.0 67.7 ± 0.4 −9.2 ± 0.2 −66.2 ± 1.8 −75.4 ± 1.8 −43.5 ± 0.9 35.8 ± 1.2 
Trypsin 2PTC_1TMY Unstable −157.7 ± 3.3 230.9 ± 3.6 73.2 ± 0.2 −8.6 ± 0.1 −61.6 ± 0.5 −70.1 ± 0.5 −39.7 ± 0.6 42.8 ± 1.1 

aUnits in kcal/mol; ΔGELE, non-solvent electrostatic potential energy; ΔGPB, electrostatic contributions to the solvation free energy calculated with the Poisson–Boltzmann equation; ΔGPBELE, electrostatic potential energy; ΔGSA, non-polar contributions to solvation free energy; ΔGVDW, van der Waals potential energy; ΔGNP, non-polar potential energy; T ΔSNM, the entropic contribution calculated with normal mode analysis to the free energy of binding; ΔGcalc, calculated binding free energy.

In the lysozyme designer complexes, the overall free energy was also found to be more favorable for the stable complex (1VFB_1TMY). The non-polar component (ΔGNP) was also more favorable in the stable complex, whereas the entropy of binding appeared to be similar for all three complexes. Interestingly, one of the unstable complexes exhibited more favorable electrostatic binding free energy (ΔGPBELE). Unlike the non-polar component, electrostatics were not found to be a strong predictor of stability.

Finally, there was also one stable complex among the three trypsin designer complexes. The overall free energy of binding reflected this observation, showing a significant more favorable binding affinity for the complex between trypsin and MBP1 (1bm8). The large difference in the overall binding affinity is reflected throughout the components of the free energy that are consistently more favorable in this complex.

Flexible backbone design

We selected four complexes from among the nine complexes (Fig. 4) that were designed with fixed backbone. We chose two stable (1VFB_1TMY and 2PTC_1BM8) and two unstable (1BRS_1BM8 and 1VFB_1TEN) complexes, and redesigned their interface using multiple conformers from MD simulations. In each case, 50 conformers collected from the MD simulation of the complex were used along with our design algorithm using Eqs (1)–(3). The final four complexes were then immersed in explicit solvent and subjected to 30 ns of MD. The r.m.s. deviation over the course of the 30 ns trajectories (Fig. 7) revealed three of the complexes remained stable, and one, namely the complex between lysozyme and tenascin (1VFB_1TEN), was unstable and appeared to dissociate.

Fig. 7

Mean r.m.s. deviation for the four complexes that were designed with flexible backbone.

Fig. 7

Mean r.m.s. deviation for the four complexes that were designed with flexible backbone.

To gain further insight into the changes that multiconformer design introduced, we performed extensive free energy calculations and the results are provided in Table VII. As expected, the unstable complex revealed significantly less favorable non-polar component of the free energy ΔGNP. Interestingly, the unstable complex had the most favorable electrostatics ΔGPBELE and configurational entropy T ΔSNM. This is counterintuitive as it would be expected that these components should not be favorable. We then compare the components of the free energy for the fixed and flexible backbone designs listed in Tables V and VII, respectively. In the case of the complexes that became stable when flexible backbone design was used, namely 1BRS_1TMY, the flexible backbone design complex exhibited a significantly more favorable overall free energy of binding from 41 to −20 kcal/mol (ΔGcalc). The non-polar component of the free energy, ΔGNP, also showed significant improvement, and likely was the key factor that contributed to the stability of this complex. More favorable electrostatics were also observed, but entropy appeared to be less favorable in the flexible backbone designer complex, suggesting a possible mechanism for improving binding. For the complex that remained unstable, a slight improvement in the overall affinity was observed. The other components also showed improvement, including the non-polar and entropy, but electrostatics were less favorable in the flexible backbone designer complex. Comparison of the r.m.s. deviation for fixed and flexible backbone designs of 1VFB_1TEN, however, shows that the flexible backbone showed greater structural stability, as the mean r.m.s. deviations remain within 4 Å over the course of the trajectories (Fig. 7), whereas the fixed backbone complexes reaches 7 Å (Fig. 5B). Finally, it is worth noting that among the two complexes that were initially stable for the fixed backbone and remained stable in the flexible backbone design (1VFB_1TMY and 2PTC_1BM8), one (1VFB_1TMY) showed significant improvement in the free energy of binding. Although the overall affinity was less favorable in the case of 2PTC_1BM8, the nearly identical non-polar component likely was the basis for the stability that was exhibited by the complex.

Table VII

Free energy calculations for complexes designed with flexible backbonea

Receptor PDB Stability ΔGELE ΔGPB ΔGPBELE ΔGSA ΔGVDW ΔGNP T ΔSNM ΔGcalc 
Barnase 1BRS_1TMY_Redesign Stable −417.8 ± 2.0 459.5 ± 1.7 41.7 ± 0.1 −11.8 ± 0.0 −96.3 ± 0.4 −108.1 ± 0.4 −46.2 ± 0.4 −20.3 ± 1.1 
Trypsin 2PTC_1BM8_Redesign Stable 77.3 ± 1.1 −20.5 ± 1.3 56.9 ± 0.1 −10.5 ± 0.0 −93.4 ± 0.3 −103.9 ± 0.3 −37.4 ± 0.5 −9.6 ± 0.8 
Lysozyme 1VFB_1TEN_Redesign Unstable −424.1 ± 2.5 463.1 ± 2.0 39.0 ± 0.1 −8.1 ± 0.0 −58.2 ± 0.3 −66.3 ± 0.4 −31.6 ± 0.3 4.3 ± 0.7 
Lysozyme 1VFB_1TMY_Redesign Stable 147.9 ± 2.1 −105.2 ± 1.8 42.7 ± 0.1 −9.3 ± 0.0 −79.5 ± 0.3 −88.8 ± 0.3 −33.9 ± 0.4 −12.2 ± 0.8 
Receptor PDB Stability ΔGELE ΔGPB ΔGPBELE ΔGSA ΔGVDW ΔGNP T ΔSNM ΔGcalc 
Barnase 1BRS_1TMY_Redesign Stable −417.8 ± 2.0 459.5 ± 1.7 41.7 ± 0.1 −11.8 ± 0.0 −96.3 ± 0.4 −108.1 ± 0.4 −46.2 ± 0.4 −20.3 ± 1.1 
Trypsin 2PTC_1BM8_Redesign Stable 77.3 ± 1.1 −20.5 ± 1.3 56.9 ± 0.1 −10.5 ± 0.0 −93.4 ± 0.3 −103.9 ± 0.3 −37.4 ± 0.5 −9.6 ± 0.8 
Lysozyme 1VFB_1TEN_Redesign Unstable −424.1 ± 2.5 463.1 ± 2.0 39.0 ± 0.1 −8.1 ± 0.0 −58.2 ± 0.3 −66.3 ± 0.4 −31.6 ± 0.3 4.3 ± 0.7 
Lysozyme 1VFB_1TMY_Redesign Stable 147.9 ± 2.1 −105.2 ± 1.8 42.7 ± 0.1 −9.3 ± 0.0 −79.5 ± 0.3 −88.8 ± 0.3 −33.9 ± 0.4 −12.2 ± 0.8 

aUnits in kcal/mol; ΔGELE, non-solvent electrostatic potential energy; ΔGPB, electrostatic contributions to the solvation free energy calculated with the Poisson–Boltzmann equation; ΔGPBELE, electrostatic potential energy; ΔGSA, non-polar contributions to solvation free energy; ΔGVDW, van der Waals potential energy; ΔGNP, non-polar potential energy; T ΔSNM, the entropic contribution calculated with normal mode analysis to the free energy of binding; ΔGcalc, calculated binding free energy.

In sum, we have implemented a flexible backbone approach toward the redesign of the interface of existing protein–protein complexes. We found that incorporating receptor flexibility led to stable complexes in all three cases. In addition, the trypsin–BPTI designer complex exhibited significantly greater affinity than the native complex due to more favorable electrostatic, non-polar and entropy components. We then grafted the surface of three proteins, namely MBP1, CheY and tenascin, to bind to barnase, trypsin and lysozyme receptors. This was first done with fixed backbone and resulted in at least one stable complex for each receptor. This was particular encouraging given the difficulty that we had previously encountered in designing ligands that bind to barnase in a stable manner. Our results show that the non-polar component of the free energy to be consistently more favorable for the stable complexes. The electrostatic component, however, in some cases was more stable for the unstable complexes, confirming our previous results that favorable electrostatic interactions may not be sufficient to predict stability. The configurational entropy was another component that did not show any particular trend. Among the nine designer complexes, we then picked two stable and two unstable and redesigned the interfaces by applying our design algorithm on multiple conformers collected from MD simulations of the complexes. We found that one of the unstable complexes became stable, whereas the other gained additional stability. In each case, the flexible design led to complexes that showed significantly greater affinity than the complexes that were designed with fixed backbone.

Supplementary data

Supplementary data (colour images for figures 2, 4 and 6) are available at PEDS online.

Funding

The research was supported by the INGEN grant from the Lilly Endowment, Inc. Computer time on the Big Red supercomputer at Indiana University is funded by the National Science Foundation under Grant No. ACI-0338618l, OCI-0451237, OCI-0535258 and OCI-0504075 as well as by Shared University Research grants from IBM, Inc. to the Indiana University. This research was supported in part by the Indiana METACyt Initiative. The Indiana METACyt Initiative of Indiana University is supported in part by Lilly Endowment, Inc.

References

Ali
M.H.
Taylor
C.M.
Grigoryan
G.
Allen
K.N.
Imperiali
B.
Keating
A.E.
Structure
 , 
2005
, vol. 
13
 (pg. 
225
-
234
)
Binz
H.K.
Amstutz
P.
Pluckthun
A.
Nat. Biotechnol.
 , 
2005
, vol. 
23
 (pg. 
1257
-
1268
)
Bonvin
A.M.J.J.
Curr. Opin. Struct. Biol.
 , 
2006
, vol. 
16
 (pg. 
194
-
200
)
Butterfoss
G.L.
Kuhlman
B.
Annu. Rev. Biophys. Biomol. Struct.
 , 
2006
, vol. 
35
 (pg. 
49
-
65
)
Chaudhury
S.
Gray
J.J.
J. Mol. Biol.
 , 
2008
, vol. 
381
 (pg. 
1068
-
1087
)
Chen
R.
Li
L.
Weng
Z.
Proteins
 , 
2003
, vol. 
52
 (pg. 
80
-
87
)
Clavel
T.
Germon
P.
Vianney
A.
Portalier
R.
Lazzaroni
J.C.
Mol. Microbiol.
 , 
1998
, vol. 
29
 (pg. 
359
-
367
)
Dall'Acqua
W.
, et al.  . 
Biochemistry
 , 
1998
, vol. 
37
 (pg. 
7981
-
7991
)
Fu
X.
Apgar
J.R.
Keating
A.E.
J. Mol. Biol.
 , 
2007
, vol. 
371
 (pg. 
1099
-
1117
)
Green
D.F.
Dennis
A.T.
Fam
P.S.
Tidor
B.
Jasanoff
A.
Biochemistry
 , 
2006
, vol. 
45
 (pg. 
12547
-
12559
)
Humphris
E.L.
Kortemme
T.
Structure
 , 
2008
, vol. 
16
 (pg. 
1777
-
1788
)
Jain
M.
Kamal
N.
Batra
S.K.
Trends Biotechnol.
 , 
2007
, vol. 
25
 (pg. 
307
-
316
)
Janin
J.
Henrick
K.
Moult
J.
Eyck
L.T.
Sternberg
M.J.
Vajda
S.
Vakser
I.
Wodak
S.J.
Proteins
 , 
2003
, vol. 
52
 (pg. 
2
-
9
)
Jiang
L.
, et al.  . 
Science
 , 
2008
, vol. 
319
 (pg. 
1387
-
1391
)
Johnson
E.S.
Annu. Rev. Biochem.
 , 
2004
, vol. 
73
 (pg. 
355
-
382
)
Keskin
Z.
Gursoy
A.
Ma
B.
Nussinov
R.
Chem. Rev.
 , 
2008
, vol. 
108
 (pg. 
1225
-
1244
)
Kortemme
T.
Baker
D.
Proc. Natl Acad. Sci. USA
 , 
2002
, vol. 
99
 (pg. 
14116
-
14121
)
Kortemme
T.
Joachimiak
L.A.
Bullock
A.N.
Schuler
A.D.
Stoddard
B.L.
Baker
D.
Nat. Struct. Mol. Biol.
 , 
2004
, vol. 
11
 (pg. 
371
-
379
)
Krol
M.
Chaleil
R.A.
Tournier
A.L.
Bates
P.A.
Proteins
 , 
2007
, vol. 
69
 (pg. 
750
-
757
)
Kuhlman
B.
Dantas
G.
Ireton
G.C.
Varani
G.
Stoddard
B.L.
Baker
D.
Science
 , 
2003
, vol. 
302
 (pg. 
1364
-
1368
)
Lensink
M.F.
Mendez
R.
Curr. Pharm. Biotechnol.
 , 
2008
, vol. 
9
 (pg. 
77
-
86
)
Liang
S.
Liu
S.
Zhang
C.
Zhou
Y.
Proteins
 , 
2007
, vol. 
69
 (pg. 
244
-
253
)
Liang
S.D.
Li
L.W.
Hsu
W.L.
Pilcher
M.N.
Uversky
V.
Zhou
Y.Q.
Dunker
A.K.
Meroueh
S.O.
Biochemistry
 , 
2009
, vol. 
48
 (pg. 
399
-
414
)
Lippow
S.M.
Tidor
B.
Curr. Opin. Biotechnol.
 , 
2007
, vol. 
18
 (pg. 
305
-
311
)
Lippow
S.M.
Wittrup
K.D.
Tidor
B.
Nat. Biotechnol.
 , 
2007
, vol. 
25
 (pg. 
1171
-
1176
)
Liu
Y.
Liu
Z.
Androphy
E.
Chen
J.
Baleja
J.D.
Biochemistry
 , 
2004
, vol. 
43
 (pg. 
7421
-
7431
)
Liu
S.
Liu
S.Y.
Zhu
X.L.
Liang
H.H.
Cao
A.N.
Chang
Z.J.
Lai
L.H.
Proc. Natl Acad. Sci. USA
 , 
2007
, vol. 
104
 (pg. 
5330
-
5335
)
Marvin
J.S.
Hellinga
H.W.
Proc. Natl Acad. Sci. USA
 , 
2001
, vol. 
98
 (pg. 
4955
-
4960
)
Massova
I.
Kollman
P.A.
Perspect. Drug Discov. Des.
 , 
2000
, vol. 
18
 (pg. 
113
-
135
)
Montclare
J.K.
Schepartz
A.
J. Am. Chem. Soc.
 , 
2003
, vol. 
125
 (pg. 
3416
-
3417
)
Ogata
K.
Jaramillo
A.
Cohen
W.
Briand
J.P.
Connan
F.
Choppin
J.
Muller
S.
Wodak
S.J.
J. Biol. Chem.
 , 
2003
, vol. 
278
 (pg. 
1281
-
1290
)
Pichler
A.
Knipscheer
P.
Oberhofer
E.
van Dijk
W.J.
Korner
R.
Olsen
J.V.
Jentsch
S.
Melchior
F.
Sixma
T.K.
Nat. Struct. Mol. Biol.
 , 
2005
, vol. 
12
 (pg. 
264
-
269
)
Ray
M.C.
Germon
P.
Vianney
A.
Portalier
R.
Lazzaroni
J.C.
J. Bacteriol.
 , 
2000
, vol. 
182
 (pg. 
821
-
824
)
Razeghifard
R.
Wallace
B.B.
Pace
R.J.
Wydrzynski
T.
Curr. Protein. Pept. Sci.
 , 
2007
, vol. 
8
 (pg. 
3
-
18
)
Rosenberg
M.
Goldblum
A.
Curr. Pharm. Des.
 , 
2006
, vol. 
12
 (pg. 
3973
-
3997
)
Schreiber
G.
Fersht
A.R.
J. Mol. Biol.
 , 
1995
, vol. 
248
 (pg. 
478
-
486
)
Sia
S.K.
Kim
P.S.
Proc. Natl Acad. Sci. USA
 , 
2003
, vol. 
100
 (pg. 
9756
-
9761
)
Smith
G.R.
Sternberg
M.J.
Bates
P.A.
J. Mol. Biol.
 , 
2005
, vol. 
347
 (pg. 
1077
-
1101
)
Srinivasan
J.
Cheatham
T.E.
III
Kollman
P.A.
Case
D.A.
J. Am. Chem. Soc.
 , 
1998
, vol. 
120
 pg. 
9401
 
Vita
C.
, et al.  . 
Proc. Natl Acad. Sci. USA
 , 
1999
, vol. 
96
 (pg. 
13091
-
13096
)
Yang
W.
Wilkins
A.L.
Ye
Y.
Liu
Z.R.
Li
S.Y.
Urbauer
J.L.
Hellinga
H.W.
Kearney
A.
van der Merwe
P.A.
Yang
J.J.
J. Am. Chem. Soc.
 , 
2005
, vol. 
127
 (pg. 
2085
-
2093
)
Yin
H.
Slusky
J.S.
Berger
B.W.
Walters
R.S.
Vilaire
G.
Litvinov
R.I.
Lear
J.D.
Caputo
G.A.
Bennett
J.S.
DeGrado
W.F.
Science
 , 
2007
, vol. 
315
 (pg. 
1817
-
1822
)
Yosef
E.
Politi
R.
Choi
M.H.
Shifman
J.M.
J. Mol. Biol.
 , 
2009
, vol. 
385
 (pg. 
1470
-
1480
)
Zondlo
N.J.
Schepartz
A.
J. Am. Chem. Soc.
 , 
1999
, vol. 
121
 (pg. 
6938
-
6939
)
Edited by Valerie Daggett