S.Liang et al. more than 0.1 A˚ 2 upon complexation. Solvent-accessible surface area is calculated as described by Zou et al.68 Side- chains of other residues and the backbones are fixed. Side- chain conformations of interface residues are optimized by Monte Carlo simulated annealing simulation as described in Refs. 65, 62. First, the rotamers for the interface residues are initialized by random selections. Second, an interface residue is randomly selected and the frequency to select the residue is proportional to the number of rotamers pos- sessed by the residue. Third, a rotamer for the residue is selected at random and the interaction energy between the rotamer and the rest of the protein E nsiedwechain is calculated using Eq. (1). The change of rotamer for the residue is accepted if the energy value is decreased (E nsiedwechain < E osliddechain). Otherwise, the change or the move is accepted with probability exp[(E soilddechain E nsiedwechain)/T]. Each cycle has 50N trials of rotamer substitutions or 5N successful substitutions, whichever comes first (N is the number of interface residues). The annealing temperature is set to 10 initially and reduced by a factor of 0.8 after each cycle. A total of 15 cycles of annealing are conducted. Energy minimization. The resulting structure from side- chain modeling is further minimized by CHARMM.14 This is to remove backbone steric clashes and discrete errors of the rotamers. Fifty steps of adopted-basis Newton Raphson minimizer (ABNR) are applied. All charged resi- dues are in their charged states and distance-dependent dielectric constant is used to calculate electrostatic energy. CHARMM 19 parameters are used throughout minimiza- tion. Here, the solvent effect is only approximated by a distance-dependent dielectric constant. In the ZDOCK decoys, only heavy atoms have coordi- nates. All polar hydrogen atoms are added to decoy sets with the program REDUCE.67 Nonpolar hydrogen atoms are ignored. To reduce the possible effect of overoptimiza- tion and computational time, only 3/4 of sidechain-opti- mized decoys with high binding affinity are selected for energy minimization. We discard 1/4 of decoys with the lowest predicted binding affinities after sidechain model- ing. This is because those decoys usually contain severe atomic clashes and could be overoptimized if they are used in energy minimization. EMPIRE score function We develop the EMPIRE score function by extending the sidechain energy described earlier for the evaluation of binding affinity. Equation (1) becomes I E ¼ 0:143 ScIontact þ 0:724 Voverlap I þ 1:72 EhIbond þ 28:6 EeIlec 0:0467 DSpIho þ 0:0042 DSpIhi I 30 I þ 1:14 DðFphiÞ þ 7:95 Vexclusion: 2 ð Þ Here, each term is evaluated between two binding proteins rather than between a given sidechain and a protein in 246 PROTEINS sidechain optimization. For example, ScIontact is the con- tacting surface area between two proteins and DSIpho (DSpIhi) is the difference between buried hydrophobic (hydrophilic) solvent accessible surface areas before and after the binding. Both sidechain and mainchain atoms are included in calculation. Unlike Eq. (1) for protein design, the terms f1,f2,Nssbond, and DGref make no contri- bution to binding affinity and are not considered. How- ever, we need to add a reference value of the protein–pro- tein interaction to Eq. (2) instead of amino acid reference value in Eq. (1) to calculate binding free energy. Similar to the role of amino acid reference value (DGref) in Eq. (1), the interface reference value is the nonspecific interaction energy between the two proteins in a refined decoy com- plex. Because the amino-acid reference values strongly cor- relate with the size of amino acid residues,62 we hypothe- size that the interface reference value is proportional to the buried solvent accessible surface area of the interface. That is, DG ¼ c1EI þ c2SI þ c3; 3 ð Þ where c1, c2, and c3 are to-be-determined constants (see Results) and SI is the buried solvent accessible surface area of interface obtained with a solvent probe of radius 1.4 A˚ . RESULTS EMPIRE score function We obtain the three constants (c1, c2, and c3) of the ref- erence state [Eq. (3), methods] by multiple linear-regres- sion analysis between EI, SI, and DGexp. The training set (1ppf, 3sgb, 1jhl, 2ptc, 1cho, 1acb, 2kai, 1vfb, 1mlc, 1nmb, 1brs, and 1ycs) are selected from 75 complexes collected by Conte et al.69 Here, a complex is selected for training if (1) its structure has a resolution of 2.5 A˚ or higher, (2) it does not undergo disorder–order transitions upon associa- tion, (3) its binding affinity70 is known, and (4) its buried interface area is between 1200 and 1600 A˚ 2. The medium size of binding interface is used here because significant conformational changes between complexed and free components are observed69 for large interfaces of 2000– 4660 A˚ 2. Thus, they are not suitable for the purpose of obtaining the reference state. We first attempted to use Eq. (3) for linear regression with the above-described training set. However, we found that DGexp does not have a significant correlation with SI. This leads to a nearly negligible c2. Thus, we employ EI as the dependent variable in the regression, instead. This yields the equation: EI ¼ 3.09 DGexp 0.0295 SI þ 3.77. Then, a simple inversion yields c1 ¼ 0.324, c2 ¼ 0.00955, and c3 ¼ 1.22. The energy unit is kilocalories per mole. The correlation coefficient between the calculated and experimental binding free energies is 0.73. The fact that c2 > 0 indicates the unfavorable contribution to association. DOI 10.1002/prot