Protein-Protein Docking with F2Dock 2.0 and GB-Rerank

Motivation Computational simulation of protein-protein docking can expedite the process of molecular modeling and drug discovery. This paper reports on our new F2 Dock protocol which improves the state of the art in initial stage rigid body exhaustive docking search, scoring and ranking by introducing improvements in the shape-complementarity and electrostatics affinity functions, a new knowledge-based interface propensity term with FFT formulation, a set of novel knowledge-based filters and finally a solvation energy (GBSA) based reranking technique. Our algorithms are based on highly efficient data structures including the dynamic packing grids and octrees which significantly speed up the computations and also provide guaranteed bounds on approximation error. Results The improved affinity functions show superior performance compared to their traditional counterparts in finding correct docking poses at higher ranks. We found that the new filters and the GBSA based reranking individually and in combination significantly improve the accuracy of docking predictions with only minor increase in computation time. We compared F2 Dock 2.0 with ZDock 3.0.2 and found improvements over it, specifically among 176 complexes in ZLab Benchmark 4.0, F2 Dock 2.0 finds a near-native solution as the top prediction for 22 complexes; where ZDock 3.0.2 does so for 13 complexes. F2 Dock 2.0 finds a near-native solution within the top 1000 predictions for 106 complexes as opposed to 104 complexes for ZDock 3.0.2. However, there are 17 and 15 complexes where F2 Dock 2.0 finds a solution but ZDock 3.0.2 does not and vice versa; which indicates that the two docking protocols can also complement each other. Availability The docking protocol has been implemented as a server with a graphical client (TexMol) which allows the user to manage multiple docking jobs, and visualize the docked poses and interfaces. Both the server and client are available for download. Server: http://www.cs.utexas.edu/~bajaj/cvc/software/f2dock.shtml. Client: http://www.cs.utexas.edu/~bajaj/cvc/software/f2dockclient.shtml.


Details on Affinity Function Computations
This section adds detail on the FFT based implementation and parameter selection of three affinity functions: shape complementarity, electrostatics and interface propensity (or hydrophobicity).

Shape Complementarity
First we present the new skin-core definitions in detail.

Skin-core Definition and Weighting
-Floating Receptor (i.e., Stationary Molecule A) Skin: Unlike the traditional approach [1] the receptor skin layer does not touch the receptor van der Waals surface (i.e., there is a gap between the skin and the core regions), and also the radius of those skin atoms differ from that in the traditional approach. The plots in Figure S2 justify this skin-core gap which, for each of 60 rigid-body test cases chosen from Zlab Benchmark Suite 2.0 [2], plots average distance and standard deviation of ligand atom centers from the receptor van der Waals surface when the ligand is bound to the receptor. We have considered only those ligand atoms that lie within 3Å from the receptor surface in bound condition. For the bound-bound case (when both molecules are crystallized together in the complex form) the average distance lies between 1.6Å and 2.0Å, while for unbound-unbound case (when the two molecules are crystallized separately) it lies between 1.4Å and 2.0Å most of the time. The standard deviation lies between 0.55Å and 0.75Å for the bound-bound case, and between 0.6Å and 1Å for the unbound-unbound case. In our experiments pseudo skin atoms of radius 1.1Å with their centers placed at a distance of 1.7Å from the receptor vdW surface (so there is a gap of 0.6Å) performed much better than the traditional approach where pseudo atoms of radius 1.4Å are placed touching the receptor vdW surface. The weights assigned to receptor skin atoms are computed based on the curvature of the skin around that atom (as opposed to the constant weights assigned to skin atoms in the traditional approach [1]). Such weighting encourages convex-concave and concave-convex interfaces as opposed to large flat interfaces.
We know that the centers of molecule A's skin atoms are placed at a constant distance d from the receptor vdW surface S. Let us enlarge S outward by distance d so that it passes through the skin atom centers, and let V be the volume enclosed by the enlarged surface. Let B k be a ball of radius δ at the center of skin atom k, and B in k = B ∩ V and B out k = V \ B in k . We define: Then the weight assigned to skin atom k of molecule A is c SC A,k = ρ(k). We approximate ρ(k) from the 3D grid embedding of molecule A along with its skin region, and δ = 3.6 produced good results in our experiments.
-Depth dependent Weighting of Receptor Core: The core atoms of molecule A are assigned weights using an increasing function of depth (distance of the atom center from the surface of A). Such weightings discourage deeper core-core overlaps more compared to shallower ones by making it difficult to compensate for deep core-core overlaps with wide skin-skin overlaps.
The core atom weights are computed based on atomic layers. Any atom which contributes to the solvent excluded surface (SES) is considered an exposed atom. The exposed atoms of a molecule lie in layer 1. When atoms in layer 1 are removed, the new exposed atoms make layer 2. In general, when atoms in layers 1 to l ≥ 1 are removed, the exposed atoms among the remaining ones make layer l + 1. The atomic layer of core atom k of A is denoted by L A (k), and the weight assigned to that atom is c SC Hence overall, the weight assigned to atom/pseudo-atom k of molecule A: -Thicker Skin of the Ligand (i.e., Moving Molecule B) with Skin-Core Overlap: Since in the traditional approach [1] the ligand skin is defined using its surface atoms, the skin thickness varies and can be too thin in some areas. It turns out that a double layer ligand skin with skin-core overlap performs better than a single skin layer disjoint from the core. The weights assigned to the atoms of B are as follows. otherwise.

FFT based formulation
Now suppose the user-specified relative weights for various types of overlaps are as follows: w ss = reward for (unit) skin-skin overlap, w cc = penalty for (unit) core-core overlap, and w sc = reward/penalty for (unit) skin-core overlap.
Let A ′ denote molecule A with its grown skin layer. If the weight assigned to an atom/pseudoatom k (of molecule P ∈ {A ′ , B} ) is c SC k = c SC k,Re + i · c SC k,Im , then the affinity function for shape complementarity is: The overall shape complementarity score for translation t and rotation r is: with T and ∆ being the translation and the rotation operator, respectively.

Electrostatics (E).
In [3] Gabb et. al. described a simplified model for electrostatics which allows efficient FFT-based computation of the term during docking search. We used this simple electrostatics model in the original version of F 2 Dock 2.0 [1]. The first protein's electric potential is computed and matched against the charges in the other. Charge assignments are made using PDB2PQR [4]).
Two affinity functions f E A and f E B are defined for molecule A and B, respectively.
where, q k is the Coulombic charge on atom k, δ(x) is the Kronecker delta function with value 1 at ||x|| = 0, and 0 everywhere else, g E A,k (x) = g E B,k (x) = 1 and E(x) is the distance dependent dielectric constant [3] as given below.

FFT based formulation
In the current version of F 2 Dock 2.0, we define for P ∈ A, B, where, γ > 0 is a constant, and blobbiness β = 2.3. In our experiments the new definition of g E P,k (x) with γ = 3.4 performed much better than the original Gabb et al. formulation. Smoothing with the Gaussian function as above has the effect of reduced discretization error on the grid.
The overall electrostatics score for translation t and rotation r is: , and w E is the user-specified weight given to electrostatics interaction.
1.3 Interface Propensity (IP) and Hydrophobicity (HP). F 2 Dock 2.0 scores the interfaces between molecules A and B using the per-residue interface propensity values computed in [5] based on the relative frequencies of different residues in the interfaces of a set of 63 protein-protein complexes from [6]. Let IP (R) denote the natural logarithm of the interface propensity value of a residue R. Alternatively, the user can choose to use per-residue Hydrophobicity values from [7]. Note that, Hydrophobic residues are more likely to be found on the interface, and hence Hydrophobicity and Interface propensity model the same chemical property and hence are used as substitutes in F 2 Dock 2.0, as opposed to using both at the same time.
The IP values for the 20 amino acid residues lie between -0.38 (for ASP) and 0.83 (for TRP). A residue with a higher IP value is likely to occur more frequently in a protein-protein interface than one with a lower IP value. Hydrophobic residues Phe, Met, Ile, Leu and Val, polar aromatic residues Trp, Tyr and His, and the charged residue Arg favor interface locations, and they have positive IP values. Additionally, Cys and Asn also have positive IP . The IP value of an atom is set to the IP value of the residue it belongs to.
Suppose A + B t,r is a docking pose obtained by rotating molecule B by r and translating by t. Let iAtom + t,r (P ) and iAtom − t,r (P ) denote the set of atoms in the interface of P ∈ {A, B} in this docking pose that have positive and negative IP values, respectively. Also let iAtom s t,r (A, B) = iAtom s t,r (A) + iAtom s t,r (B), for s ∈ {+, −}. Then we assign the following interface propensity score to the docking pose: where IP ǫ = max IP (R)<0 IP (R).

FFT based formulation
We show below how we approximate IP -score t,r (A, B for A and B, respectively, and We set the radius of each atom of Q to its van der Waals radius, while the van der Waals radius of each atom of Q ′ is extended by a non-negative constant r ext 1 . We assign weight c IP,Q P,k to each atom k ∈ P = {Q, Q ′ } = {A, B} as follows: Then we define: and 0 otherwise. The overall interface propensity score for translation t and rotation r (i.e., approximated IP -score t,r (A, B)) is: where, w IP is the user-specified weight given to interface propensity.

Details on Filters
In this section we describe the implementations of various filters available in F 2 Dock 2.0. Each filter has a preprocessing step (e.g., construction of octrees and initialization of other relevant data structures) which is executed only once before F 2 Dock 2.0 enters the rotations loop in phase I, and so not included in the running times of the filters which are executed in each iteration of the loop. Also octrees constructed for the initial configuration of molecule B can be easily transformed on-the-fly so that it correctly matches any rotated configuration of B by applying the corresponding rotation matrix on the node coordinates (e.g., node center) and atoms/points stored in them.

Lennard-Jones Filter.
We approximate the Lennard-Jones (LJ) potential between molecules A and B t,r given by the following expression.
where r ij is the distance between atoms i ∈ A and j ∈ B t,r , constants a ij and b ij depend on the type (e.g., C, H, O, etc.) of the two atoms involved. For any fixed pair of atom types a ij and b ij are fixed, and are calculated from the Amber force field using well depths µ and equivalence contact distances of homogeneous pairs r eqm as follows (assuming X = atomT ype(i ∈ A) and Y = atomT ype(j ∈ B t,r )) [8,9].
We assume that X, Y ∈ {C, H, N, O, P, S}.
Our approach is to identify potential docking solutions based on a cutoff potential. Assuming that the docking solution can not have steric clashes at the interface, we can filter (or penalize) all poses with positive LJ potential. However, though this assumption is not unreasonable for boundbound docking, in the case of bound-unbound and unbound-unbound docking we need to allow some steric clashes. Hence we soften the potential by reducing the r eqm,XX values by a constant factor δ LJ which has the effect of reducing the inter-atomic clash distances 2 by the same factor. As a result docking poses with negative LJ values can now have some clashes, and the maximum number of allowed clashes will depend on the value of δ LJ . In our experiments δ LJ = 0.3 produced good results. Observe is the number of atoms in molecule A (resp. B). Instead we use our octree-based fast multipole-type algorithm described in [10] to obtain a 1 + ǫ approximation of the potential in is a user-defined approximation parameter.

Interface Area Filter.
We sample weighted quadrature/integration points from the surface of each molecule as described in [11,10]. Sum of the weights of all such integration points is equal to the molecular surface area up to quadrature approximation error, and each integration point can be treated as the center of a small surface patch with its weight being the patch area. Each such patch π with center c π and area a π is represented by the tuple c π , a π .
We then evaluate the following expression as an approximation of the interface area between molecules A and B t,r . where, and µ IA is a user-defined constant.
Using our algorithm described in [10] based on octrees [12] and our Dynamic Packing Grid (DPG) data structure [13]

Interface Propensity Filter.
As in Section 2.2 we assume that we have a patch decomposition of the surface of each molecule P ∈ {A, B}. However, each patch π is now represented as a triple c π , a π , h π where the additional value h π is a measure of interface propensity of π in the following sense. We assume that each atom k ∈ P is assigned a interface propensity value h k as in Section 1.3, and h π = 1 |S P,π | k∈S P,π h k , where S P,π = {k|(k ∈ P ) ∧ (dist(c k , c π ) < νr k )}, c k is the center of atom k, r k is its van der Waals radius, and ν is a user-defined constant.
This filter computes IP V ec(A, B t,r ) which is defined as the following vector: where IP P 1 (P 2 ) is as defined in Section 2.2.
This vector can be computed simultaneously with the approximate interface area in Section 2.2 incurring negligible overhead. The algorithm is given in [10].
Once we have IP V ec(A, B t,r ) = v 1 , v 2 , v 3 , v 4 , we approximate IP -score t,r (A, B) defined in Section 1.3 as follows: . We also observed that rewarding docking poses with larger positive interface propensity weighted area v 2 + v 4 generally improves the ranks of near native poses. However, rewarding based on the product of IP -score t,r (A, B) and v 2 + v 4 seems to work even better in practice. This composite term which we define as IP -product t,r (A, B) below, rewards a docking pose with higher v 2 + v 4 value more than one with lower such value when both have the same IP -score.
The interface propensity filter in F 2 Dock 2.0 penalizes a docking pose provided its IP -score is below a user-defined threshold µ IP . The IP -product value, after weighting it by a user-specified weight w ′ IP , is added to the overall score of the pose.

Clash Filter.
Two atoms a ∈ A and b ∈ B with van der Waals radii r a and r b , respectively, are said to be in a clash provided the distance between their centers is smaller than α(r a +r b ), where α is a user-defined positive constant. F 2 Dock 2.0 counts the total number of atomic clashes between molecules A and B t,r as follows.
Clash t,r (A, B) = |{b|(b ∈ B t,r ) ∧ ∃ a∈A (r a,b < α(r a + r b ))}| Direct computation of this number requires O (M A M B ) time. So we use an octree-based spatial subdivision approach that is similar to the interface area approximation algorithm used in Section 2.2 and has similar running time.
F 2 Dock 2.0 penalizes a docking pose A + B t,r provided Clash t,r (A, B) > µ C , where µ C is a user-defined constant. A value of around 10 for µ C seems to work well in practice.   [14] based on user preference with Table  III being the default. The contact preferences were derived as follows. Two residues are considered to be in contact if the distance between their C β atoms (C α for Gly) is less than 6 Å. The normalized number of contacts between residue types i and j is defined as Q ij = C ij / k<l C kl , where C ij is the actual number of contacts observed between residue types i and j. The expected number of contacts W ij is obtained by assuming that there are no contact preferences between residues of different types. Then the likelihood of contact between residues types i and j is defined as (Table  III of [14]): G ij = 10 log Table IV of [14]) is obtained by replacing Q ij with:

Residue-Residue Contact Filter.
where V i is the volume of residue i. Given a docking pose A + B t,r , we identify all residue-residue contacts at the interface of the two molecules using a fast algorithm similar to the one used in Section 2.3, and compute the sum of all positive and negative G ij values (or G ij (v) values if chosen by the user) denoted by G + and G − , respectively. Then we compute the following ratio: if the user chooses volume-normalized contact preferences). F 2 Dock 2.0 penalizes a docking pose if its RC-score is below a user-specified constant µ R C. Without volume normalization µ R C = 3.0 seems to work well in practice. Docking poses can also be penalized based on user-specified lower and upper bounds for G + and G − , respectively.

Glycine Filter.
We mark the oligopeptides with the properties described above on the enzyme surface, and for any given docking pose count the number of these motifs occurring at the interface using a fast algorithm similar to the one used for counting atomic clashes in Section 2.4.

Antibody-Antigen Contact Filter.
We say that an atom a is in close neighborhood of another atom b provided the distance between the centers of a and b is at most 2(r a + r b ). Given a potential antibody-antigen docking pose, F 2 Dock 2.0 computes three quantities: N L1∪H1 , N L3 and N H3 , denoting the number of antigen atoms that are in the close neighborhood of any atom in the antibody regions CDR-L1/CDR-H1, CDR-L3 and CDR-H3, respectively. These counts are obtained using a fast algorithm similar to the one used in Section 2.4. F 2 Dock 2.0 penalizes a docking pose provided N L1∪H1 < µ L1∪H1 or N L3 < µ L3 or N H3 < µ H3 , otherwise it adds w AC (N L1∪H1 + N L3 + N H3 ) to the total score, where µ L1∪H1 , µ L3 , µ H3 and w AC are user-specified quantities. F 2 Dock 2.0 uses µ L1∪H1 = 100, µ L3 = µ H3 = 150 and w AC = 10 by default. The CDR (Complementarity Determining Region) loops are identified using the method described in [15].

Overall Cost of Filtering.
Summing up the cost of all filters described above for a single run, we obtain Assuming that each filter is applied on at most N F configurations, the total time taken by all filters is where N R is the number of samples in the rotations space, the running time reduces to

Solvation Energy Based Reranking
The solvation energy E sol consists of the energy to form cavity in the solvent (E cav ), the solutesolvent van der Waals interaction energy (E vdw(s-s) ), and the electrostatic potential energy change due to the solvation (also known as the polarization energy, E pol ) [16,17,18,19].
The first two terms are often modeled as [16,20] where p is the solvent pressure, V is the molecular volume, A i is the solvent accessible surface area of atom i and γ i is its solvation parameter, ρ 0 is the bulk density, and u (att) i is the van der Waals dispersive component of the interaction between atom i and the solvent.
The last term, E pol , can be approximated using Generalized Born (GB) theory [21].
where τ = 1 − 1 ǫ , and R i is the effective Born radius of atom i. GB-rerank approximates each of these terms as described in the following sections, and reranks the list of top docking poses produced by F 2 Dock 2.0 based on the resulting ∆E sol values. In order to approximate ∆E sol , GB-rerank precomputes the E sol values for molecules A and B, and then computes E sol for each docking pose.

Approximating E pol (Polarization Energy).
We first estimate all Born radii of the molecule/complex using our fast approximation scheme described in [10] followed by another of our algorithms from [10] for estimating E pol from the approximated Born radii. Given a molecule/complex consisting of M atoms, both algorithms run in O 1 ǫ 3 M log M time using O (M ) space, where ǫ > 0 is an approximation parameter that provides a speed-accuracy tradeoff. The smaller the value of ǫ the more accurate the algorithms are, and the larger the ǫ value the faster they run.
GB-rerank uses either the r 4 -or the r 6 -approximation of R i (see [10]) based on user preference. By default, it uses the r 6 -approximation [22] which shows better accuracy for spherical solutes as well as for proteins [23].
where Γ is the molecular surface, n(r) is the unit outward normal on the molecular surface at r, and x i is the center of atom i.
We obtain the following discrete surface formulation of Equation (2) by applying the divergence theorem and Gaussian quadrature.
where, the r k 's are N = O (M ) Gauss quadrature/integration points [11] on the molecular surface, and w k is a weight assigned to r k in order to achieve higher order of accuracy for small N . In order to avoid the slowdown due to the repeated computation of quadrature points for each docking pose A + B t,r , F 2 Dock 2.0 precomputes the set of quadrature points Q A and Q B for A and B, respectively, and chooses all quadrature points from Q A ∪ Q B t,r that are not on the interface of A and B t,r for use during the Born radii estimation of A + B t,r .

Approximating E vdw(s-s) (Dispersion Energy).
The solute-solvent van der Waals interaction energy (also known as dispersion energy) is modeled as [16,20]: where ρ 0 is the bulk density, and u (att) i is the van der Waals dispersive component of the interaction between atom i ∈ [1, M ] and the solvent which is given as follows.
If R i is the Born radius of atom i calculated using the r 6 -approximation (i.e., Equation 2/3), then Equation 4 can be rewritten as: Therefore, E vdw(s-s) can be approximated in O (M ) time once the Born radii of all atoms are available.

Approximating E cav (Cavity Forming Energy).
Instead of computing E cav for A, B and each A + B t,r separately, we approximate ∆E cav with the buried surface area of A + B t,r which is approximated using the same algorithm used for interface area filter in Section 2.

Overall Cost of Reranking.
Assuming that GB-rerank is applied on N G docking poses, its total running time is where N R is the number of samples in the rotations space, and so the running time reduces to

Additional Results on the Effects of Various Affinity Functions and Filters
Figures S3 to S10 track the changes in the rank of the top hit, RMSD of the top hit, number of hits in top 1000, and the minimum RMSD in the top 1000 predictions, respectively, as the various F 2 Dock 2.0 options are activated one after another on the rigid-body test cases from Zlab benchmark 2.0 [2]. When we activate Lennard-Jones filter, clash filter and proximity clustering after shape complementarity we get hits for 4 new test cases, and the rank of the top hit improves for 15 more (see Figure S3(top)). However, we also lose hits in top 1000 for 3 test cases, and the rank of the top hit degrades for one test case. Overall, the application of these filters and clustering seem largely beneficial. The best results are obtained for enzyme-inhibitor/enzyme-substrate complexes, as for more than 50% of these complexes rank of the top hit improve.
The impact on the RMSD of the top hit is not clear though it has improved for 5 test cases in the enzyme-inhibitor/enzyme-substrate group (see Figure S5(top)). Similar is the impact on the minimum RMSD in the top 1000 predictions (see Figure S9(top)).
The number of hits in the top 1000 predictions is reduced in majority of the cases (see Figure  S7(top)) as the two filters penalize predictions that are close to the target complex but have steric clashes. The proximity clustering step also penalizes poses that are very similar to poses with higher scores.
When electrostatics is turned on we get hits in top 1000 for 9 test cases for which we did not have a single hit before, and for 14 other cases rank of the top hit improve (see Figure S3(middle)). However, we lose all hits for 1 test case, and for 4 others rank of the top hit degrades.
Though there does not seem to be any particular trend in the change of the RMSD of the top hit except for the 9 cases with new first hits as mentioned above (see Figure S5(middle)), for antibodyantigen and antigen-bound antibody complexes the minimum RMSD in the top 1000 seems to have generally improved (see Figure S9(middle)). For 22 test cases the number of hits in the top 1000 predictions has improved while for 5 cases this number has dropped (see Figure S7(middle)).
The FFT-based interface propensity scoring is activated next which improves the rank of the top hit for 30 test cases (i.e., for around 50% of all cases) among which 7 cases did not have a single hit before (see Figure S3(bottom)). Among these 7 cases with new first hits 5 are antibodyantigen or antigen-bound antibody complexes, and none are enzyme-inhibitor or enzyme-substrate. Ranks have degraded for 5 test cases. Minimum RMSD in top 1000 has improved for 20 test cases degrading for only 1 (see Figure S9(bottom)), and RMSD of top hit has also degraded for only 1 test case (see Figure S5(bottom)). The number of hits has increased for 28 test cases and decreased for 7 (see Figure S7(bottom)).
The interface propensity filter is turned on next. It improves the rank of the top hit for 25 complexes, and degrades for 5 (see Figure S4(top)). For 3 test cases we did not have a single hit in top 1000 before among which 2 are antibody-antigens. The number of hits has increased for 16 cases, and decreased for 5 (see Figure S8(top)). The impacts on the RMSD of the top hit and the minimum RMSD in top 1000 are not significant in general (see Figures S6(top) and S10(top)).
The residue-residue contact filter which is activated next improves the rank of the top hit for 27 test cases, and degrades for none (see Figure S4(middle)). The enzyme-inhibitor and enzymesubstrate complexes seem to have benefited the least from this filter. The impact on the number of hits in top 1000 is not significant and is mixed (see Figure S8(middle)). The RMSD of the top hit and the minimum RMSD in top 1000 do not change except for the 3 cases for which we get new first hits (see Figures S6(middle) and S10(middle)).
Next we apply the antibody contact filer and the Glycine filter. The antibody contact filter improves the rank of the top hit for 9 antibody-antigen and antigen-bound antibody test cases, and degrades for 3, while the Glycine filter slightly improves the same for 4 enzyme-inhibitor/enzymesubstrate complexes (see Figures S4(bottom)). The antibody contact filter improves the number of hits in the top 1000 predictions for 5 complexes (see Figure S8(bottom)).
All other changes are insignificant.