Affinity- and Specificity-Enhancing Mutations Are Frequent in Multispecific Interactions between TIMP2 and MMPs

Multispecific proteins play a major role in controlling various functions such as signaling, regulation of transcription/translation, and immune response. Hence, a thorough understanding of the atomic-level principles governing multispecific interactions is important not only for the advancement of basic science but also for applied research such as drug design. Here, we study evolution of an exemplary multispecific protein, a Tissue Inhibitor of Matrix Metalloproteinases 2 (TIMP2) that binds with comparable affinities to more than twenty-six members of the Matrix Metalloproteinase (MMP) and the related ADAMs families. We postulate that due to its multispecific nature, TIMP2 is not optimized to bind to any individual MMP type, but rather embodies a compromise required for interactions with all MMPs. To explore this hypothesis, we perform computational saturation mutagenesis of the TIMP2 binding interface and predict changes in free energy of binding to eight MMP targets. Computational results reveal the non-optimality of the TIMP2 binding interface for all studied proteins, identifying many affinity-enhancing mutations at multiple positions. Several TIMP2 point mutants predicted to enhance binding affinity and/or binding specificity towards MMP14 were selected for experimental verification. Experimental results show high abundance of affinity-enhancing mutations in TIMP2, with some point mutations producing more than ten-fold improvement in affinity to MMP14. Our computational and experimental results collaboratively demonstrate that the TIMP2 sequence lies far from the fitness maximum when interacting with its target enzymes. This non-optimality of the binding interface and high potential for improvement might characterize all proteins evolved for binding to multiple targets.


Introduction
Virtually all functions in the cell are regulated through cascades of protein-protein interactions (PPIs). Some biological processes cause activation of several parallel PPI pathways that frequently intertwine with each other. At the crossroads of such pathways lie proteins that are capable of interacting with a number of different partners and hence are called multispecific proteins [1]. Due to their central role in PPI networks, multispecific proteins are crucial to cell survival and their malfunction inevitably leads to disease. Thus, unraveling the atomic-based principles for binding multispecificity is not only interesting for basic biology but also valuable for the studies directed at finding new therapeutics that target various PPIs. Binding interface sequences of multispecific proteins are under evolutionary pressure to provide favorable interactions for various partners that in some cases share little sequence and structure homology. These sequences are a compromise required for accommodating multiple targets and thus cannot be optimal for interactions with each individual target protein. We postulate that binding interface sequences of multispecific proteins lie far from the fitness maximum for each individual interaction and thus could be further improved through mutations. In other words, mutations that enhance binding affinity should be frequent in multispecific PPIs. Moreover, such mutations are likely to narrow down binding specificity of multispecific proteins towards a particular target or a set of targets.
To test this hypothesis, we chose a representative multispecific protein, Tissue Inhibitor of Metalloproteinases 2 (TIMP2). TIMP2 is one of four similar proteins in humans (TIMP1, 2, 3 and 4) that regulate a family of more than twenty-six homologous enzymes, Matrix Metalloproteinases (MMPs) and the related ADAMs families [2][3][4]. MMPs play a major role in degradation of the extracellular matrix and participate in many important biological processes such as embryonic development, organ morphogenesis, bone remodeling and others. On the other hand, imbalance in MMP activity is associated with a diverse set of diseases including arthritis, cardiovascular diseases, neurological disorders, fibrosis, and cancer [5]. MMPs are multi-domain proteins that differ in domain architecture and substrate preferences [6] but all share a catalytic domain with a nearly identical active site containing a Zn 2+ ion. High-resolution structures have been solved for a number of MMPs alone and in complex with TIMPs [7][8][9][10][11][12][13]. [Murphy, 2011 #643] These structures reveal that TIMPs bind directly to the catalytic zinc ion at the active site of the enzyme, shielding it from the solvent. The interaction is conveyed mostly through the TIMP N-terminal domain (N-TIMP) consisting of ,125 residues. N-TIMP is a potent inhibitor of various MMPs and thus has been repeatedly used as a substitute for the full-length protein in various experimental studies [14]. N-TIMP binds to MMPs mostly through four contiguous regions ( Figure 1A). The first region includes six N-terminal residues that come in close proximity to the enzyme active site and coordinate the catalytic Zn 2+ through the N-terminal Cys. Besides the Nterminal region, three additional N-TIMP loops (35-42, 66-72, and 97-99 in N-TIMP2) participate in direct interactions with MMPs ( Figure 1A).
MMPs are synthesized in an inactive form. They could be activated by other MMPs and inactivated upon binding of TIMPs. Each of the four known TIMPs binds most of the MMPs with slightly different affinities, ranging from 10 211 -10 29 M. In addition, some TIMPs, such as TIMP-2, can participate in the activation path of certain MMPs, through binding to another MMP type [15]. TIMP/MMP interactions hence present a complicated regulatory network with connections that are not fully understood. Rational manipulation of this network through mutations could help to elucidate precise functional roles of various MMPs and facilitate development of selective inhibitors for each MMP type. TIMPs present an attractive scaffold for design of such inhibitors and hence have been a subject of various mutational studies. Previous studies demonstrated that certain substitutions at positions 2, 4, and 68 of TIMP2 can strongly affect its relative affinity for different MMPs [16][17][18]. In another study, a single mutation T98L that stabilizes TIMP1 in the bound conformation was shown to produce an impressive specificity shift towards MMP-14 relative to other MMP types [19,20]. More recently, phage display technology was used to probe a large number of possible mutations in the N-TIMP2 binding interface and to engineer a variant that binds to MMP1 with a nanomolar affinity while losing its affinity to MMP3 and MMP14 [21].
In contrast to previous studies, our goal was to obtain a more comprehensive picture of TIMP/MMP interactions and to locate positions on TIMP where affinity-and specificity-enhancing mutations could be identified with high probability. For this purpose, we first generated computational binding landscapes of N-TIMP2/MMP14 and N-TIMP2/MMP9 interactions by predicting effects of all single mutations in the N-TIMP2 binding interface on its affinity to these two enzymes. We validated some of our predictions experimentally by constructing a number of N-TIMP2 mutants and measuring their affinity to these two enzymes. We extended our computational studies to six additional MMPs for which structural models for interactions with TIMP2 could be constructed. Both computational and experimental results point to the suboptimal nature of the N-TIMP2 binding interface sequence and possibility of affinity and specificity improvement through various mutations. These results are in agreement with our hypothesis that multispecific proteins are not optimized for a particular binding partner and could be reengineered to be more selective interaction partners and inhibitors.

Materials and Methods
Model construction for N-TIMP2/MMP complexes We created models for structures of MMP/TIMP2 complexes for those MMPs that have their X-ray structure available only in the unbound form (MMP1 (PBD 3SHI), MMP2 (PDB 1RTG), MMP3 (PDB 1B3D), MMP7 (PDB 1MMQ), and MMP9 (PDB 1L6J)). For this purpose, we first superimposed the unbound structure of a particular MMP on the structure of the MMP14/N-TIMP2 complex (PDB 1BUV). We next removed from the structure all MMP residues that do not belong to the catalytic domain. This initial superimposed structure of the MMP/N-TIMP2 complex was then refined using the RosettaDock server [22]. The best output structure from the RosettaDock server was used as an input for the saturation mutagenesis protocol.

Computational saturation mutagenesis
An In silico saturation mutagenesis protocol was applied on the N-TIMP2 binding interface using the structure of the N-TIMP2/ MMP-14 (PDB code 1BUV) complex [8], the N-TIMP2/MMP13 (PDB code 2E2D) complex [9], and the N-TIMP2/MMP10 complex (PDB code 4ILW) and the models constructed for the N-TIMP2/MMP complex. Only the N-terminal TIMP2 domain was used in the calculations (residues 1-127 of TIMP2). The metal ions Ca 2+ and Zn 2+ were not considered in the calculations. The N-TIMP2 binding interface was defined as all the residues that are within 4 Å from the MMP in the N-TIMP2/MMP14 structure and included residues 1-4, 6, 14, 35, 36, 38, 40, 42, 66, 68, 69, 70, 71, 97, 99, and 100-101. From this set we excluded positions that are very close to the catalytic zinc ion (positions 1-3, and 100-101) and positions that coordinate a calcium ion (position 36) since effects on catalysis and interactions with ions could not be adequately modeled by our protocol. The remaining fourteen residues were scanned using the saturation mutagenesis protocol described in Sharabi et al [23]. Briefly, for each of the fourteen positions, we performed 18 calculations where the considered position on N-TIMP2 was either kept WT or was replaced with another amino acid, all except for Pro, Cys, and Gly. During the calculation, the interface residues as well as residues in direct contact with the interface were repacked and the energy of the N-TIMP2/MMP complex was calculated for the WT and for the mutated complex. We then separated the two chains and calculated the energy of each single chain. The intermolecular energy was calculated by subtracting the energies of the single chains (N-TIMP2 and MMP) from the total energy of the complex. DDG bind was calculated by subtracting the intermolecular energy of the WT complex from the intermolecular energy of the mutant. Finally, the obtained DDG bind was normalized according to a linear equation obtained in our previous work where correlation between various experimental and computed DDG bind values was tested [24]. Rotamer libraries used for design were based on the backbone dependent library of Dunbrack and Karplus [25] with additional rotamers expanded by one standard deviation around their mean x 1 and x 2 values. For the calculations, we used ORBIT software with the energy function optimized by our group for design of protein-protein interactions [24]. The energy function contained terms that describe Van der Waals attractive and repulsive interactions, hydrogen bond interactions, electrostatic interactions, and surface-area-based solvation (see [24] for the exact description of the energy function,). The lowest-energy rotameric conformation of each mutant was found using the Dead-End Elimination theorem [26,27]. Finally, we color-coded each mutation according to its DDG bind value in kcal/mol: DDG bind $1.5 -red; 0.5#DDG bind ,1.5 -yellow; 20.5#DDG bind ,0.5 -green and DDG bind #20.5 -blue. Mutations that were predicted to destabilize unbound N-TIMP2 or an unbound MMP by more than 2 kcal/mol were considered potentially deleterious for N-TIMP2 folding and were colored in gray if predicted to improve DDG bind .

Evaluating N-TIMP2 position tolerance and specificity potential
We evaluated tolerance of each N-TIMP2 binding interface position for mutations based on the results of the saturation mutagenesis protocol for DDG bind prediction. For this purpose, we replaced each color saturated mutagenesis figure a score: 21 for blue, 0 for green, 1 for yellow, and 2 for red mutations. Gray mutations were not incorporated in the calculation. We calculated the average score over all mutations at a single binding interface position for one MMP and assigned positions into three classes according to the score: Score #0.2 R tolerant, 0.2,Score#1 R semi-tolerant, Score .1 R non-tolerant. To evaluate the potential of a particular mutation to narrow down binding specificity, we compared DDG bind predictions for one particular mutation among the eight MMP types. For each particular mutation, we calculated the average score and its standard deviation over all MMPs. A mutation with standard deviation greater than 1 was considered beneficial for enhancing binding specificity over all eight MMPs and was marked by a star. In addition, we calculated an average score and standard deviation over all mutations for each N-TIMP2 position.

MMP enzymes
Catalytic domains of MMP14 and MMP9 were expressed recombinantly and purified as published before [12].

Expression and refolding of the N-TIMP2 mutants
Genes for the N-TIMP-2 mutants were generated by the Transfer PCR protocol [28] starting from the plasmid pET-28atimp-2-HISX6 containing the gene for the WT N-TIMP2 (residues  All experiments were repeated 3 times and average and standard deviation for K d is reported 2 The specificity shift was defined as the fold of affinity enhancement towards MMP14, divided by the fold of affinity enhancement towards MMP9. 3 The specificity is considered predicted correctly if the calculated difference between DDG bind for the N- Binding affinity determinations using the enzyme activity essay The synthetic fluorogenic MMP substrate MCA-Pro-Leu-Gly-Leu-Dpa-Ala-Arg-NH 2 ?TFA [where MCA is (7-methoxycoumarin-4-yl)acetyl; Dpa is N-3-(2,4-dinitrophenyl)-L-2,3-diaminopropionyl; and TFA is trifluoroacetic acid] was purchased from GenScript Inc., (Piscataway, NJ) and used to assay enzyme activity. Samples with various concentrations of the N-TIMP2 mutant were pre-incubated with either MMP9 (at 0.

Results
Mapping computational binding landscapes for N-TIMP2/MMP14 and N-TIMP2/MMP9 interactions Each protein-protein interaction can be characterized by a binding landscape that represents changes in protein-protein binding affinity due to point mutations. To generate the binding landscape of the N-TIMP2/MMP14 and N-TIMP2/MMP9 interactions, we used the computational saturation mutagenesis protocol developed in our lab [29]. This protocol scans each PPI binding interface position with all amino acids, repacks the surrounding side chains and determines the change in free energy of binding due to mutations (DDG bind ) (see Methods). As an input for the protocol, we used an X-ray structure of the N-TIMP2/ MMP14 complex and a structural model of the N-TIMP2/MMP9 complex generated from the structure of unbound MMP9 (see Methods). While usage of a structural model instead of an actual structure is bound to introduce some inaccuracies in our calculations, we were optimistic in the case of MMP9 since this enzyme exhibits high structural homology to MMP14 in the N-TIMP2 binding region with a Ca RMSD of 0.66 Å for interfacial residues ( Figure 1B). Next, we computationally scanned fourteen N-TIMP2 positions with seventeen amino acids, producing binding energy landscapes for the N-TIMP2/MMP14 and the N-TIMP2/MMP9 interactions. We did not consider mutations to Gly and Pro since those mutations are likely to result in backbone conformational changes that were not modeled by our protocol. Mutations to Cys were not considered due to their possible To better visualize the binding energy landscape of the N-TIMP2/MMP14 complex, each mutation was assigned to one of four classes according to the predicted DDG bind value and colored in blue, green, yellow, and red for affinity-enhancing, neutral, destabilizing, and highly destabilizing mutations, respectively ( Figure 2). In addition, we classified each N-TIMP2 binding interface position according to its ability to accept various mutations into tolerant, semi-tolerant and non-tolerant (see Methods). Figure 2 shows that the N-TIMP2 interface is not particularly optimized for binding to either MMP since the landscape is dominated by neutral and affinity-enhancing mutations represented by green and blue circles. For example, for MMP-14, out of fourteen considered N-TIMP2 binding interface positions, ten showed possibility of significant DDG bind improvement with mutation to at least one amino acid (Blue circles, Figure 2). For

Experimental testing of computational predictions
To determine how well our computational binding landscapes reflect the reality of the TIMP2/MMP binding energetics, we decided to validate some of the predictions experimentally. The number of the tested N-TIMP2 mutants was limited by a relatively tedious procedure for their construction that requires refolding after expression in E. Coli [14]. We hence selected thirteen N-TIMP2 single mutants, focusing on mutations that 1) were predicted to enhance binding affinity to MMP14 and 2) were predicted to enhance binding specificity towards MMP14 relative to MMP9 (Table 1). To measure binding between the N-TIMP2 mutants and MMP14/MMP9, we utilized an enzyme activity assay described previously ( Figure 3A) [30]. This assay is based on detecting fluorescence resulting from cleavage of a fluorogenic MMP substrate. High sensitivity of the assay allows us to measure binding affinities as low as 10 211 M.
Using the above assay, we determined K d s for interactions between N-TIMP2 WT and MMP14 and MMP9 to be 4.5 and 0.9 nM respectively, similar to previously published results [21]. These K d s became a point of reference for calculating DDG bind for the selected N-TIMP2 mutants. Eight mutations that were predicted to significantly improve binding affinity of N-TIMP2 to MMP14 (S4R, S4Q, V6R, I35K, N38Q, S68W/Y, H97R) proved to be affinity-enhancing experimentally (Table 1 and Figure 3B). In addition, two mutations, S4A and T99Y that were predicted as neutral or slightly destabilizing also showed improved binding affinity for MMP14. Among the affinity-enhancing mutations two, I35K and H97R, exhibited a 12-and 14-fold improvement in binding affinity towards MMP14, an impressive affinity shift for single mutations. Five of the N-TIMP2 mutants with increased affinity towards MMP14 (S4R, V6R, I35K, S68W, and H97R) also exhibited affinity enhancement towards MMP9 ( Figure 3B and Table 1). This demonstrates that affinityenhancing mutations at the N-TIMP2 binding interface could be easily found through our computational protocol. Twelve out of thirteen explored N-TIMP2 mutants produced a detectible shift in binding specificity towards MMP14 relative to MMP9 ( Table 1). Three of the mutations produced a substantial (7-9 fold) shift in binding specificity towards MMP14 relative to MMP9, including mutations I35K and H97R that also improved affinity towards this enzyme and S4E that was slightly destabilizing for the complex with MMP14 and highly destabilizing for the complex with MMP9 (Table 1). We hence conclude that both affinity-and specificity-enhancing mutations are quite frequent at the N-TIMP2 binding interface. Good consensus of our predictions with experimental results for MMP9, where no actual X-ray structure of the complex was available, gave us confidence that fairly realistic binding landscapes could be constructed with our computational saturation mutagenesis protocol using structural models of PPIs as starting points. We next tested whether our findings on low-optimality of the N-TIMP2 binding interface could be extended to additional MMP types.

Computational binding landscapes of N-TIMP2 interacting with additional MMPs
We aimed to explore N-TIMP2 computational binding landscapes for as many MMPs as possible. The limiting factor here was the structural information on the N-TIMP2/MMP complexes. Presently, three high-resolution structures of the TIMP2/MMP complexes are available (for MMP10, MMP13 and MMP14). Nevertheless, we were able to generate additional structural models of the N-TIMP2 complexes for those MMPs that have their structure solved in the unbound form including: MMP1, MMP2, MMP3, and MMP7. Together with MMP14 and MMP9, we thus analyzed computational binding landscapes for eight MMP family members. Figure 4

Discussion
Low optimality of the N-TIMP2/MMP binding landscapes Our computational and experimental findings point to the relatively low optimality of the N-TIMP2/MMP interfaces. The in silico saturation mutagenesis protocol predicted that affinity enhancement could be produced through at least three different mutations at eight different N-TIMP2 positions when interacting with MMP14 and at six different positions when interacting with MMP9 ( Figure 2). Experimentally, with only a few trials, we identified affinity-enhancing mutations at seven and five positions for MMP14 and MMP9, respectively (Table 1). Computational binding landscapes of the remaining MMPs also point to high potential for affinity improvement (Figure 4). The low optimality of the N-TIMP2/MMP interfaces is not surprising since N-TIMP2 binds to all MMP members with similar affinities and hence cannot provide favorable intermolecular interactions for each MMP. Our results are in agreement with previous computational studies of multispecific proteins whose binding interface sequence was found to be optimal for simultaneous interactions with different targets yet sub-optimal for interaction with each target on its own [31,32]. An ability to greatly improve binding affinity and specificity through only a few mutations was recently experimentally demonstrated in ubiquitin, a protein whose function is to bind to multiple targets with low affinity [33]. In contrast, our recent study on high-affinity enzyme-inhibitor complexes revealed highly optimized binding landscapes with only a handful of mutations that further increase affinity [34]. All of the above findings suggest that low optimality of the binding interface might be a general property of multispecific interactions that distinguishes them from PPIs with narrow binding specificity.

Analysis of affinity-enhancing mutations
Eight out of ten experimentally identified affinity-enhancing mutations were correctly predicted for the N-TIMP2/MMP14 interaction and four out of five mutations were correctly predicted for the N-TIMP2/MMP9 interaction, demonstrating the potential of our in silico saturation mutagenesis approach in identifying affinity-enhancing mutation and its applicability not only to crystal structures but also to structural models. Among the identified mutations two, I35K and H97R, exhibited affinity improvement of more than ten-fold, which is higher than usually observed for single mutations [35]. Both of our best affinity-enhancing mutations are substitutions to positively charged residues. This is not surprising since the MMP binding interface is slightly negatively charged ( Figure 5). In addition, both substitutions occur at positions, where no significant interaction with MMP14 occurs in the wild-type N-TIMP2/MMP14 complex while favorable intermolecular interactions are created upon substitution. Substitution of His to Arg at position 97 forms additional Van der Waals interactions and creates new intermolecular hydrogen bonds and electrostatic interactions with Asp 193 on MMP14 and with the backbone carbonyl ( Figure 6A). Mutation of an Ile at position 35 to Lys also improves interface packing and creates favorable electrostatic interactions with Asp 188 and Asp 212 on MMP14 ( Figure 6B). Other identified affinity-enhancing mutations (S68W, S68Y and T99Y) are mutations from small to aromatic amino acids that fill in the gaps in the non-optimal interface and bury additional hydrophobic area ( Figure 6C). Burial of larger hydrophobic area has been proposed as a strategy for selecting affinity-enhancing mutations in a previous study [36].

Analysis of specificity-enhancing mutations
All but one tested N-TIMP2 mutants exhibited a shift in binding specificity towards MMP14 relative to MMP9, in agreement with most of our predictions (Table 1). While the observed specificity shift was modest, combining such mutations might result in an N-TIMP2 mutant that shows many-fold preference for one MMP type. To better understand how binding specificity is conveyed at the molecular level, we analyzed all mutations in the structural context, by looking at modeled complexes between N-TIMP2 mutants and MMP14/MMP9. Our analysis showed that most such mutations created specific interactions [37,38], such as hydrogen bonds, salt bridges, or pi-pi stacking interactions with MMP14, but were unable to form similar interactions with MMP9 due to the absence of an appropriate amino acid on the enzyme side ( Table 2). For example, Tyr 99 is predicted to form a pi-pi stacking interaction with F198 on MMP14, but lacks an aromatic interaction partner on MMP9. Similarly, Gln 38 is predicted to form a hydrogen bond with Q208 on MMP14; this residue is replaced by a Gly on MMP9, thus disallowing any hydrogen bond interaction.
Based on the computational binding landscapes of N-TIMP2/ MMP interactions generated in this work, we further propose a strategy for selecting specificity-determining positions, or positions where mutations have a high potential for narrowing down binding specificity. Such positions, (e. g. positions 4, 35, 66, 69, 70, and 71) display high standard deviation in DDG bind predictions across the whole MMP family (Table S1). Interestingly, half of these positions (4,66,71) are also in contact with positions on MMP that exhibit the highest sequence variability over eight MMP types ( Figure 5). These specificity-determining positions should be the focus of experiments that rely on selection of N-TIMP mutants with narrowed specificity from large combinatorial libraries of mutants.
Furthermore, using computational binding landscapes we can predict specific mutations that narrow down N-TIMP2 binding specificity for certain MMP types (Figure 4, indicated by stars). For example, mutation V6Y is predicted to significantly destabilize N-TIMP2 interactions with MMP3, MMP9, MMP13 while at the same time stabilizing its interactions with MMP1, MMP2, MMP7. Mutation V71R on the other hand is predicted to improve interactions with MMP2, while destabilizing complexes with MMP7, MMP9, MMP13, and MMP14. Note that predictions of specific mutations are more sensitive to inaccuracies in our computational protocol compared to predictions at the position and the interface level. The latter predictions reflect the global picture and depend only slightly on the results for each individual mutation.
In summary, we demonstrated that the N-TIMP2 binding interface is not optimal for binding to various MMPs, revealing a large number of mutations that improve binding affinity towards a particular MMP type. This sub-optimality might be a general property of mutispecific PPIs that have evolved to provide reasonable affinity for a large set of targets. It is thus relatively easy to enhance binding affinity of a multispecific protein towards one particular target, and the affinity increase is frequently coupled to an increase in binding specificity.

Supporting Information
Table S1 Position Specificity Potential. Average score and standard deviation over all mutations at a single binding interface position for all MMPs is presented. (DOCX)