Polar Desolvation and Position 226 of Pancreatic and Neutrophil Elastases Are Crucial to their Affinity for the Kunitz-Type Inhibitors ShPI-1 and ShPI-1/K13L

The Kunitz-type protease inhibitor ShPI-1 inhibits human neutrophil elastase (HNE, K i = 2.35·10−8 M) but does not interact with the porcine pancreatic elastase (PPE); whereas its P1 site variant, ShPI-1/K13L, inhibits both HNE and PPE (K i = 1.3·10−9 M, and K i = 1.2·10−8 M, respectively). By employing a combination of molecular modeling tools, e.g., structural alignment, molecular dynamics simulations and Molecular Mechanics Generalized-Born/Poisson-Boltzmann Surface Area free energy calculations, we showed that D226 of HNE plays a critical role in the interaction of this enzyme with ShPI-1 through the formation of a strong salt bridge and hydrogen bonds with K13 at the inhibitor’s P1 site, which compensate the unfavorable polar-desolvation penalty of the latter residue. Conversely, T226 of PPE is unable to establish strong interactions with K13, thereby precluding the insertion of K13 side-chain into the S1 subsite of this enzyme. An alternative conformation of K13 site-chain placed at the entrance of the S1 subsite of PPE, similar to that observed in the crystal structure of ShPI-1 in complex with chymotrypsin (PDB: 3T62), is also unfavorable due to the lack of stabilizing pair-wise interactions. In addition, our results suggest that the higher affinity of ShPI-1/K13L for both elastases mainly arises from the lower polar-desolvation penalty of L13 compared to that of K13, and not from stronger pair-wise interactions of the former residue with those of each enzyme. These results provide insights into the PPE and HNE inhibition and may contribute to the design of more potent and/or specific inhibitors toward one of these proteases.


Introduction
Elastases constitute a group of serine proteases (SPs) considered as attractive therapeutic targets due to their involvement in different pathologic processes. For example, pancreatic elastase is associated with pancreatitis, whereas proteinase 3 and HNE (UNIPROTKB: P08246) are elastases, molecular dynamics (MD) simulations and Molecular Mechanics Generalized-Born (Poisson-Boltzmann) Surface Area (MM-GB(PB)SA) free energy calculations [24][25][26][27][28][29][30]. These predictions showed that the interaction, i.e, salt bridge and hydrogen bonds, between K13 of ShPI-1 and D226 of HNE (Fig 1), is crucial to the complex formation. Accordingly, its abrogation through the in silico mutation D226A precludes the binding of ShPI-1 to the mutated enzyme. We also proposed that the presence of a Thr residue at the equivalent position of PPE (T226), which does not establish strong interactions with K13, greatly disfavors the binding of the wild-type inhibitor to this enzyme (Fig 1). Furthermore, it was predicted that the stronger interaction of both elastases with ShPI-1/K13L is largely caused by the lower polar-desolvation penalty of the L13 compared to that of K13, and not by stronger pair-wise interactions of the The interaction of both elastases with ShPI-1 and ShPI-1/K13L is represented as a matrix-like scheme. The inability of ShPI-1 to interact with PPE is represented by a red cross on a hypothetical complex structure. Experimentally-determined K i values for the other complexes are shown together with their respective 3D structures in cartoon representation. The crystal structure of PPE in complex with ShPI-1/K13L complex (PDB: 3UOU) was used as a template to model the structures of the remaining complexes (see Materials and Methods below). The HNE structure was extracted from the PDB 2Z7F. It is noteworthy that a hypothetical structure of the non-existing PPE:ShPI-1 complex was also generated to predict the underlying structural and energetic factors preventing its formation in solution. The P1 site residues (K13 and L13) and the residues at position 226 of both elastases (T226 and D226) are shown in stick representation.
former with the residues at the S1 subsite of each enzyme. Overall, these results point out the importance of considering both pair-wise interactions and desolvation effects when analyzing the relative affinities of protein-protein complexes.

Materials and Methods
Prediction of protease-inhibitor 3D structures through structural alignment and in silico mutations Previous lines of evidence have demonstrated that the same binding mode is preserved in protein-protein complexes formed by close homologues with 30-40% sequence identity [25]. Therefore, the 3D structure of the PPE:ShPI-1/K13L complex (PDB:3UOU) was used here as a template to generate 3D models of HNE (~39% sequence identity with PPE, S1 Fig) in complex with both inhibitor variants, i.e., HNE:ShPI-1 and HNE:ShPI-1/K13L. Likewise, the 3D model of the complex between the wild-type inhibitor ShPI-1 and PPE was calculated in order to perform structural and energetic analyses, although previous functional studies have not detected the formation of this complex [7]. First, the 3D model of the HNE:ShPI-1/K13L complex was obtained by superimposing the HNE 3D structure extracted from the crystal structure of this protease in complex with the C-terminal domain of the secretory leukocyte protease inhibitor (PDB:2Z7F) onto the PPE chain of the template (structural alignment) using the program Modeller v9.5 (Fig 2) [31,32]. Then, the 3D models of the PPE:ShPI-1 and HNE:ShPI-1 complexes were generated by performing the in silico K13L mutation with the mutagenesis tool of Pymol v1.7.0.0 (Fig 2) [33]. K13 rotamers were selected by visual inspection. This procedure was also employed for generating various alanine point mutations at the complex interfaces, which are required for computational alanine scanning (CAS) [28].

Preparation of starting structures for energy minimization and molecular dynamics simulations
The protonation states of ionizable residues and His tautomers in the protease-inhibitor complexes were determined at pH = 7.4 with the program PDB2PQR (http://nbcr-222.ucsd.edu/ pdb2pqr_1.8/), which uses PROPKA for the prediction of pK a values [34]. All acid (Glu and Asp) and basic (Lys and Arg) residues were predicted in their respective ionic forms. His residues were always neuter and appeared either in the HID or in the HIE tautomeric form. HID tautomers were specifically predicted at positions 57, 200 and 210 of PPE, and 57 and 210 of HNE in all complexes; the remaining His residues of both enzymes being predicted as HIE tautomers. The tautomeric form of H57, a residue of the catalytic triad of SPs, was particularly checked, since the information for its protonation state is available [35]. The unique His residue of ShPI-1 and ShPI-1/K13L, i.e., H47, was predicted as HID tautomer, in agreement with the structure of ShPI-1 (PDB:1SHP) determined by Nuclear Magnetic Resonance [19]. The remaining steps necessary for energy minimization (EM) were performed with GROMACS v4.5.5 package [36]. Briefly, hydrogen atoms were added to the starting structures using the protonation states of ionizable residues predicted before. Disulfide (S-S) bonds between all Cys residue pairs were built in the correct topology. Then, a dodecahedral solvation box with edges spanning at least 1 nm from the solute surface was created around each complex. TIP3P water molecules were subsequently added and periodic boundary conditions (PBC) were settled in the limits of the solvation box. Electroneutrality was guaranteed by adding Na + and Cl − ions into the unit cells at the appropriate ratio to reach final NaCl concentrations of 0.50 mol/L and 0.05 mol/L, which are similar to those used during the experimental determination of K i values against HNE and PPE, respectively [7].

Energy minimization and molecular dynamics simulations
The protocol employed here to perform MD simulations involves prior EM and positionrestrained equilibration, as outlined by Lindahl for lysozyme in water [37]. AMBER99SB forcefield [38] was used for the calculation of forces during both EM and MD simulations, which were carried out with the mdrun program of GROMACS v4.5.5 [36]. All systems were subjected to 15000 steps of steepest descents minimization with an integration step of 0.1 nm. The maximum tolerance was set to 1000 kJ mol -1 nm -1 . Cutoff radii of 1.4 nm and 1.0 nm were established for the calculation of van der Waals and short-range electrostatic interactions, respectively. The Particle Mesh Ewald algorithm was used to handle long-range electrostatic interactions [39]. Interatomic distances were left unconstrained in all systems during EM.
Subsequently, water molecules were relaxed around the complexes by 200 ps of positionrestrained equilibration. Harmonic restraints with force constants of 1000 kJ mol -1 nm -2 were applied to all heavy atoms of the proteins. The treatment of van der Waals and electrostatic Fig 2. Workflow for the Generation of the 3D Models of HNE:ShPI-1/K13L, HNE:ShPI-1 and PPE:ShPI-1 Complexes. The structural alignment of PPE and HNE was performed with Modeller v9.5 [31,32]. The in silico L13K point mutation on the inhibitor was carried out with Pymol v.1.7.0.0 [33].
interactions was identical to that of EM. The Newton's equation of motion was solved using the leap-frog algorithm, with and integration step of 2 fs [40]. Velocity rescaling [41] and Berendsen weak coupling [42] algorithms were used to keep temperature (T) and pressure (p) constant at 298 K and 1 atm, respectively. Interatomic distances were constrained by the Linear Constraints Solver algorithm [43] and random initial velocities obeying the Maxwell-Boltzmann distribution at 298 K were assigned to each atom prior to the MD simulations.
Finally, Langevin dynamics simulations [44] were carried out for each system during 25 ns at T = 298 K and p = 1 atm. The simulation time of a hypothetical non-existing PPE:ShPI-1 complex was particularly extended up to 125 ns to increase the probability of sampling the eventual disruption of the complex interface. The Parrinello-Rahman coupling algorithm [45,46] was used to keep pressure constant and the friction coefficient (ξ) was set to 0.5 ps -1 in all systems, as recommended elsewhere [36]. The treatment of non-bonded interactions and constraints, as well as the integration step were identical to those used during the positionrestrained MD simulations. Snapshots were saved at 10 ps intervals.

Contact analysis at the complex interfaces
To obtain the interatomic contacts at the interfaces of the four complexes, their representative structures were calculated from the productive (frames collected after the equilibration time (t eq )) MD simulations by using clustering analysis with g_cluster (GROMACS v4.5.5). Van der Waals contacts between residues belonging to different protein chains were then defined using a cutoff distance of 4 Å. Additionally, hydrogen bonds were calculated using the g_hbond program (GROMACS v4.5.5), based on the following geometrical criteria: i) a distance 3.5 Å between the donor and the acceptor and ii) an acceptor-donor-hydrogen angle 30 . The time stability of hydrogen bonds was also assessed during MD simulations. Finally the formation of salt bridges between oppositely-charged residue pairs interacting within a distance of 4 Å at least in one snapshot was assessed with the salt bridge extension of Visual Molecular Dynamics v1.9.1 (VMD) [47]. The average distance between the oppositely-charged groups of the interacting residues was determined from the productive MD simulations and was used as a measure of salt bridge strength.
Assessing the stability of the PPE:ShPI-1 and HNE:ShPI-1 complexes to a small interface disruption As an alternative approach to study the differential interaction of ShPI-1 with both elastases, we assessed the stability of the starting 3D models of the PPE:ShPI-1 and HNE:ShPI-1 complexes to small disruptions of their respective interfaces. This task was accomplished by first moving the ShPI-1 molecule 3.2 Å away from the PPE molecule using VMD v1.9.1. After the perturbation, the amine group of K13 (P1) of ShPI-1 was placed at the entrance of the S1 subsite of PPE. Subsequently, the disrupted HNE:ShPI-1 complex was generated by superimposing the HNE structure onto that of PPE; thus, the same perturbed coordinates of ShPI-1 were used in both complexes. The time evolution of the disrupted systems was assessed by performing 60 ns MD simulations, following the previously-described steps and conditions. The structural changes during the MD simulations were determined by performing RMSD analysis and distance time profiles. that of the effective free energy (ΔG eff ), which usually suffices for comparing the relative affinities of a series of similar ligands for a given receptor [30,48]. The term ΔE gas , which includes the internal (ΔE int ), the van der Waals (ΔE vw ) and the electrostatic (ΔE el ) energies, is derived from the force-field equations. The term ΔG solv is further decomposed into two components, i.e., the polar-solvation free energy (ΔG GB/PB ) and the non-polar solvation free energy (ΔG SA ) [30,48]. The former is calculated through different Generalized-Born (GB) or Poisson-Boltzmann (PB) implicit-solvation models, whereas the latter is obtained by the equation: in which ΔSA(X) represents the solvent-accessible surface area variation of the solute molecule X upon complex formation, while γ and β are empiric constants whose values for GB models are almost always 0.0072 kcalÁÅ -2 Ámol -1 and 0, respectively [48,49]. Finally, TÁΔS is frequently computed by normal-mode analysis and is, therefore, the most computationally-demanding step of the MM-GB(PB)SA method [48,50]. The MMPBSA.py program of Amber12 package was used for MM-GB(PB)SA free energy calculations [48,50] after converting the GROMACS trajectories into the Amber format by using VMD v1.9.1 as described elsewhere [51]. In all cases we followed the single trajectory (ST) approach, in which the trajectories for the free enzyme and the free ligand are extracted from that of the complex [48,50]. GB OBC1 , GB OBC2 and GBn2 implicit-solvation models (igb = 2, igb = 5 and igb = 8, respectively) as well as the PB model were employed to calculate the ΔG GB/PB value of each complex [29,48,[52][53][54]. Topologies were obtained with tleap using mbondi2 and mbondi3 radii for GB OBCs and GBn2 models, respectively [48]. Salt concentrations of 0.50 mol/L and 0.05 mol/L were set for HNE and PPE complexes, respectively, and default solvent and solute dielectric constants (ε w = 78.3 and ε in = 1, respectively) and rgbmax cutoff (rgbmax = 25 Å) values were used in all GB calculations [48]. In turn, default solvent probe radius (1.4 Å), dielectric constants and grid parameters were employed to solve the PB equation [48]. Additionally, three different sets of atomic radii, i.e., mbondi2, mbondi3 and Tan and Luo pre-computed values, were utilized to calculate the solute cavity [48,50]. In all cases, the ΔSA values were determined using the Linear Combination of Pair-wise Overlaps algorithm [55]. ΔG SA was then estimated through Eq 1 by setting the values of γ and β to 0.0072 and 0, respectively [50]. Finally, ΔE gas was estimated from AMBER99SB parameters [38]. Mean values of the energy terms were obtained by averaging over the snapshots extracted every 10 ps from each productive MD simulation. In turn, t eq was estimated by the analysis of accumulated mean values of ΔG eff and root mean square deviations (RMSD) during MD simulations [56]. Snapshots were considered as statistically independent from each other, since previous works have determined a typical correlation time for MM-GB(PB)SA energy terms of~5 ps or less [57,58].
Energetically-relevant residues, i.e., warm-and hot-spots, at the interfaces of the studied complexes were predicted by using the per-residue effective free energy decomposition (prE-FED) protocol implemented in MMPBSA.py [48,50]. Of note, warm-and hot-spot residues were defined as those with a side-chain energy contribution (ΔG sc ) to the total ΔG eff value ranging from -1.0 to -0.4 kcal/mol and -1.0 kcal/mol, respectively, as defined elsewhere [59]. The per-residue free energy contribution (ΔG res ) under the ST approach is calculated as follows [28,60]: The first term in the right-hand side of Eq 2 is the sum of one half of the pair-wise van der Waals (DE ij vw ), electrostatic (DE ij el ) and polar-solvation (DG ij GB=PB ) interaction energies between atoms i and j belonging to residue res and to any other residue, respectively. The second term stands for the sum of the self-interaction energies of all atoms belonging to residue res (DG ii GB=PB ). Finally, the third and fourth terms represent the per-residue non-polar solvation free energy (ΔG SA(res) ) and entropy (TΔS res ) contributions, respectively. ΔG SA(res) is calculated through the ICOSA algorithm [28,61]; whereas TΔS res is neglected by default in the MMPBSA. py free energy decomposition protocol [50]. Therefore, a per-residue effective free energy rather than per-residue binding free energy decomposition was obtained.
Another approach used here to assess per-residue free energy contributions was the CAS protocol [48,50,60]. Briefly, Ala single-point mutations were generated at specific positions as described before and the topologies of the mutated complexes were obtained using tleap of Amber12. Subsequently, relative free energy values (ΔΔG) between the native and mutated complexes were determined using MMPBSA.py. These calculations were performed under the ST approach in which the trajectory of the mutated complex is generated from that of the native complex by simply truncating the side-chain of the residue of interest and replacing the Cγ atom by a hydrogen [50,60]. The linear correlation between the per-residue energy contributions predicted through both protocols, i.e., CAS and prEFED, was assessed by estimating the Pearson coefficient (r p ) with Mathematica v7.0 [62]. Likewise, the ability of both approaches to similarly rank the per-residue energy contributions was evaluated with the Spearman ranking coefficient (r s ).
Finally, the pair-wise effective free energy decomposition (pwEFED) protocol of MMPBSA. py was employed to calculate interaction energies between pair of residues (ΔG r1, r2 ) [28,48,50]. The ΔG r1, r2 values were calculated through the following equation [28]: Summation in Eq 3 is carried out over all atoms i and j of residues r1 and r2, respectively. Note also that the ½ factor introduced elsewhere to avoid double counting [28] has been dropped from Eq 3 in order to obtain full pair-wise interaction energies.

Calculation of polar-desolvation and screened electrostatic energies
Some MM-GB(PB)SA energy components are associated with ideal processes which do not take place during the complex formation in solution. For example, ΔE el is referred to the variation of electrostatic energy in vacuum; whereas ΔG GB/PB quantifies the energy variation due to solute-solvent electrostatic interactions arisen from the transfer of individual reactant and product molecules from the vacuum to the solvent. Therefore, we used here a modification of the traditional polar-energy components (ΔE el and ΔG GB ) of the MM-GBSA method, first proposed by Zou et al., to assess the polar-desolvation penalty of the solute molecules (ΔG ds ) and the screened electrostatic energy variation (DG sc el ) upon complex formation [63]. These energy components were calculated as follows (see S2 Fig for a derivation of ΔG ds based on a thermodynamic cycle): DG sc el ¼ It is worth saying that Eq 5 is valid only under the ST approach, since we assumed that X i;j DE ij el ¼ 0 for i and j belonging to the same molecule (either R or L). The ΔG ds and DG sc el values for the studied complexes were determined through the pwEFED protocol implemented in the MMPBSA.py program [48,50]. Briefly, the output of this protocol was processed to retrieve the self (i, i)-and cross (i, j)-energy values for all residues in each complex, which were then added to obtain the total values of ΔG ds and DG sc el through Eqs 4 and 5, respectively.

Structural analysis of the protein-inhibitor interfaces during the MD simulations
Analysis of the interfaces of the three existing complexes. The 3D models of HNE in complex with each inhibitor variant, ShPI-1 and ShPI-1/K13L, were generated by using the crystal structure of the PPE:ShPI-1/K13L complex as a suitable template (Fig 2) and were subsequently used as the starting structures for MD simulations. The stability of the three existing complexes during the MD simulations was monitored by calculating the instantaneous RMSD values for different atom sets, which were relatively stable and modularly small (<2.5 Å) after 5 ns (= t eq ) in all cases (Fig 3). Of note, stable RMSD time profiles were obtained for the interface heavy atoms and, particularly, for those of the inhibitor's P1 site (Fig 3), thereby suggesting the good complementarity of the complex interfaces during the simulation time.
To validate the in silico structural analysis carried out here, the van der Waals contacts ( 4 Å) and hydrogen bonds occurring at the interface of the crystal structure of PPE in complex with ShPI-1/K13L were first compared to those at the corresponding representative structure obtained from the productive (t!t eq ) MD simulation (S1 Table). As shown in the table, most of the protease-inhibitor interactions observed in the crystal structure are also found in the representative structure of the complex. The main difference involved Q192 of PPE which underwent a side-chain rearrangement during the MD simulation, thereby forming a hydrogen bond with C12 at the P2 site of ShPI-1/K13L (S3 Fig). The overall correspondence between both structures was also inferred from the small global RMSD value (1.44 Å) calculated for their respective heavy atoms.
According to the structural analysis of the three existing complexes, most of the proteaseinhibitor interactions, i.e. van der Waals contacts (S2 Table), hydrogen bonds (S3 Table) and salt bridges (S4 Table), comprised residues at positions P3-P3' of the inhibitor's primary binding loop and the complementary S3-S3' subsites of both elastases. Remarkably, the residues at P3, P1, P2' and P3' sites of both inhibitor variants establish close van der Waals contacts with significantly-variable subsites of HNE and PPE in terms of their residue composition (S2 Table). To a lesser extent, the inhibitor's secondary binding loop (positions P19'-P24') is also involved in the complex formation (S2, S3 and S4 Tables).
As previously observed in similar complexes of canonical inhibitors and SPs [10,13,14], the interactions at the S1:P1 interfaces are predominant within the three protease-inhibitor interfaces analyzed here (S2 and S3 Tables). These interfaces showed both conserved and differential polar-interaction patterns among them (Fig 4). For example, the hydrogen bonds S195(N): P1(O) and G193(N):P1(O), typical of the substrate-like binding mechanism [64], were detected with high occupancies in the three complexes (Fig 4 and S3 Table). The hydrogen bonds with S214 and H57 were also found in all the interfaces, but their occupancies were rather low and variable among the complexes (S3 Table), suggesting a subordinate role of this interaction in the complex formation. The most striking difference involving the S1:P1 interfaces was detected in the HNE:ShPI-1 complex, in which two hydrogen bonds and a salt bridge are formed between D226 of HNE and K13 of ShPI-1, while similar interactions do not occur in the other two complexes (Fig 4 and S3 and S4 Tables). Hydrogen bonds also occur outside the S1:P1 interfaces of the studied complexes in agreement with previous reports for inhibitors following a substrate-like mechanism [13] (S1 Text).
Interestingly, our results also suggest that residue E44(P31') outside the inhibitor's primary and secondary binding loops forms a hydrogen bond and a salt bridge with R36 of HNE (S3 and S4 Tables and S5 Fig). However, in PPE the nearest positively-charged residue (R61) lies farther from E44; thereby precluding the formation of both the hydrogen bond and the salt bridge (S4 Table and  Analysis of a hypothetical PPE:ShPI-1 interface. The starting structure of a hypothetical PPE:ShPI-1 complex was subjected to a 125 ns MD simulation. The time profiles of RMSD values calculated for distinct sets of protein heavy atoms indicate the occurrence of conformational changes at the complex interface during the MD simulation ( Fig 5). As revealed by visual inspection of different frames, these conformational changes comprised the complex interface reorganization leading to the further exit of the K13 side-chain from the S1 subsite. In fact, the side-chain of this residue sampled two main conformations during the MD simulation i) an 'in' conformation inside the S1 subsite of PPE and ii) an 'up' conformation at the entrance of this subsite. It is worth noting that the nomenclature PPE:ShPI-1in and PPE:ShPI-1up will be used hereinafter to refer to the whole complexes bearing the 'in' and 'up' conformations of Understanding the Affinity of PPE and HNE for ShPI-1 and ShPI-1/K13L K13 side-chain, respectively. As shown in Fig 5 the 'in conformation' of K13 transiently exited the S1 subsite until the 'up' conformation was reached at t%62 ns, but no further evidence of instability of the complex interface was observed up to 125 ns. Notwithstanding, the occurrence of complete interface disruption at longer simulation times cannot be discarded.
The analysis of interactions similar to that performed for the three existing complexes was carried out for a hypothetical PPE:ShPI-1 complex. The PPE:ShPI-1in and PPE:ShPI-1up complexes were analyzed independently from two different time intervals of the MD simulation, each sampling one of the K13 side-chain conformations (from 5 ns to 25 ns, and from 70 ns to Understanding the Affinity of PPE and HNE for ShPI-1 and ShPI-1/K13L 90 ns, respectively). The van der Waals contacts across the PPE:ShPI-1in interface were similar to those at the PPE:ShPI-1/K13L interface (S2 and S5 Tables). However, the PPE:ShPI-1up complex displayed a remarkable reduction of van der Waals contacts at the Sn:Pn (n = 1, 2, 3, 4 and 5) interfaces, whereas no noticeable differences were observed at the Sn':Pn' interfaces of both conformations (S5 Table).
The hydrogen bond and salt bridge patterns at the PPE:ShPI-1in interface are similar to those of the PPE:ShPI-1/K13L complex (S3, S4 and S6 Tables). Remarkably, the K13 side-chain did not establish stable polar interactions within the S1 subsite of PPE (Fig 6 and S6 Table). Indeed, T226 of PPE only forms a very weak hydrogen bond with the ε-amine group of K13, whereas D226 of HNE forms stable hydrogen bonds and a salt bridge with the latter group, as previously stated (Figs 4 and 6 and S3, S4 and S6 Tables). On the other hand, the PPE:ShPI-1up interface possesses a reduced number of hydrogen bonds compared to that of the PPE: ShPI-1in complex, especially, at the S3:P3 and S2:P2 interfaces. It is worth noting that the hydrogen bond S195(N):K13(O), conserved throughout the other interfaces analyzed here, and characteristic of the substrate-like interaction mechanism, is abrogated in the 'up' conformation ( Fig 6B). Overall, a decrease in both the van der Waals and polar interactions at the complex interface, mainly involving the Pn side of the inhibitor binding loop and the complementary enzyme subsites, was observed upon the transition from the 'in' to the 'up' conformation. In addition, the ε-amine group of K13 is not stabilized through polar interactions with the PPE residues in any conformation. conformation of the PPE:ShPI-1 complex in which the side-chain of K13 lies within the S1 subsite of PPE was sampled (see structures at t = 0 and t = 20 ns). From t%43 ns to t%62 ns, a reorganization of the complex interface occurred in order to provide space for K13 exit (see structure at t = 60 ns). Finally, at t%62 ns the side-chain of K13 changed to an 'up' conformation and remained in that conformation up to 125 ns (see structure at t = 125 ns). Different outcomes upon small disruption of the interfaces of ShPI-1 in complex with PPE and HNE. The MD simulation of the PPE:ShPI-1 complex showed the exit of the K13 side-chain from the S1 subsite of the enzyme; however, an eventual dissociation of the complex was not sample up to 125 ns (Fig 5). Instead of conducting longer simulations, we decided to assess whether ShPI-1 remains bound to the enzyme upon an artificial disruption of the complex interface. The same exact protocol was carried out for the HNE:ShPI-1 complex, which is expected to remain associated upon slight disruption (positive control).
The outcome of this experiment is presented in Fig 7. As can be seen, ShPI-1 readily dissociates from PPE during the 60 ns MD simulation of the disrupted complex. This suggests, in turn, that the association of ShPI-1 to PPE is an unfavorable process, since a nearly-effective collision did not lead to complex formation. In contrast, ShPI-1 remained bound to HNE upon disruption (Fig 7), thereby suggesting the stability of this complex. Hence, ShPI-1 shows a differential ability to interact with PPE and HNE, in agreement with the experimental data [7]. In addition, the distance between the T226(CG2)/D226(OD1) atom of PPE/HNE and the K13(NZ) atom of ShPI-1 as a function of time was calculated (Fig 7). As can be observed, the distance between K13(NZ) and D226(OD1) decreases during the MD simulation until it reaches a plateau at~2.7 Å, indicating the propensity of these oppositely-charged groups to restore the sat bridge interaction after the interface disruption. Conversely, the distance between T226(CG2) and K13(NZ) increases during the MD simulation, consistently with their inability to establish strong pair-wise interactions (S6 Table). Hence, we propose that the formation of the PPE:ShPI-1 complex may be impaired by an unfavorable association process, disregarding the possibility of a slow dissociation rate of the complex from hypothetical bound states. Moreover, the salt bridge between D226 of HNE and K13 of ShPI-1 may be essential to the complex formation, whereas T226 of PPE seems to be irrelevant for ShPI-1 binding. Understanding the Affinity of PPE and HNE for ShPI-1 and ShPI-1/K13L Effective free energy calculations and correlation with experimental data The large difference in the affinity of ShPI-1 and ShPI-1/K13L for PPE cannot be straightforwardly deduced from the differential interaction patterns at the interfaces of both complexes (Figs 4C and 6, and S2, S3, S4, S5 and S6 Tables). Furthermore, the HNE:ShPI-1 complex has more favorable overall interactions than the other three (Fig 4A and 4B, S2, S3 and S4 Tables), but its K i value (2.35Á10 −8 M) is higher than those of the HNE:ShPI-1/K13L (1.3Á10 −9 M) and the PPE:ShPI-1/K13L complex (1.2Á10 −8 M) [7]. These discrepancies demonstrate that, in principle, it is not possible to establish a qualitative relation between the number and/or nature of the interface contacts and the binding affinity of the complexes, as has been pointed out elsewhere [65,66]. Therefore, free energy calculations are needed to determine the energetic factors underlying the relative affinity values of these systems.
According to the instantaneous ΔG eff values calculated during the MD simulations of the three existing complexes, relatively stable accumulated mean values of ΔG eff were reached after the first 5 ns (Fig 8A, 8B and 8C). Since a similar behavior was observed in the RMSD time series (Fig 3), we considered 5 ns as a suitable t eq value. Therefore, the ΔG eff mean values were calculated using the last 20 ns of each productive MD simulation. Note, however, that an abrupt instability in the time profile of the PPE:ShPI-1 complex was observed at t%62 ns, which corresponds to the exit of K13 side-chain from the S1 subsite (compare Figs 5 and 8D). Hence, for this particular complex, the free energy calculations were carried out independently for the same two MD simulation intervals employed before in the structural analysis.
The ΔG eff mean values for each complex, estimated through different GB and PB models, are shown in S7 Table; whereas the calculated relative free energy values (ΔΔG calc ) and K i ratios for all complexes are shown in Table 1. In particular, the comparison of the calculated (ΔΔG calc ) and experimental (ΔΔG exp ) relative free energy values for the interaction of both inhibitor variants with a given elastase, revealed that all the implicit-solvation models consistently predicted the higher affinity of ShPI-1/K13L for PPE, whereas only the GB OBC2 model predicted the higher affinity of this inhibitor variant for HNE (Table 1A). In fact, the ΔΔG calc value for the HNE complexes obtained with this GB model is only 0.92 kcal/mol lower than the experimental value, which is an acceptable estimation error, especially when considering that the relative entropy contribution (TÁΔΔS) was excluded from the ΔΔG prediction. According to the ΔΔG calc value for the PPE complexes obtained through the GB OBC2 model (Table 1A), a K i value of 2.18Á10 8 M -1 for the binding of ShPI-1 to PPE is obtained from the experimental K i value of the PPE:ShPI-1/K13L complex and the ΔΔG calc value. Hence, based on the previous K i value, there must be not complex formation under the experimental conditions employed in the inhibition assays, which agrees with the reported data [7].
Additionally, the ΔΔG calc values and K i ratios for the binding of each inhibitor variant to both elastases were estimated using the corresponding ΔG eff values (Table 1B). As can be observed from the table, all implicit-solvation models correctly predicted the higher affinity of ShPI-1 for HNE. However, none of the GB models was able to correctly rank the affinities of ShPI-1/K13L for both elastases (Table 1B). Apart from the inaccuracies of the implicitsolvation models, the previous result may be influenced by the exclusion of the TÁΔΔS term Understanding the Affinity of PPE and HNE for ShPI-1 and ShPI-1/K13L from the ΔΔG calc calculations, an invalid assumption for systems bearing large structural differences [48]. In fact, in contrast to the ΔΔG calc values shown in Table 1A, the systems compared in Table 1B comprise two distinct enzymes, which, in turn, possess several point variations in the residue composition of their binding sites (see S1 Fig for a comparison of the S1 subsites of PPE and HNE). In spite of that, the ΔΔG calc value for ShPI-1/K13L in For ΔΔG calculations an average ΔG eff value of PPE:ShPI-1 complex obtained from those of the two conformations (S7 Table) was employed. f The GB model used for ΔΔG calculation is indicated by the value of the igb variable between parentheses. g The PB model used for for ΔΔG calculation is indicated between parentheses, where pb1, pb2 and pb3 stand for the PB model using Tan and Luo, mbondi2 and mbondi3 atomic radii, respectively. h No PPE:ShPI-1 complex formation has been experimentally detected, hence, ΔΔG exp value is expected to be a large and negative number, and the K i ratio must tend to zero. Understanding the Affinity of PPE and HNE for ShPI-1 and ShPI-1/K13L complex with each enzyme predicted by the GB OBC2 model is again in better agreement with the corresponding ΔΔG exp value than those obtained with the GB OBC1 , GBn2 and the three PB variants (Table 1B). Interestingly, a TÁΔΔS value of 5.30 kcal/mol would explain the difference between the ΔΔG calc value obtained with the GB OBC2 model and the ΔΔG exp value for the previous complexes. In turn, the same TÁΔΔS value would not qualitatively affect the prediction of the relative affinity of ShPI-1 for both enzymes, i.e., a large relative affinity would be obtained, ΔΔG calc = -22.70 kcal/mol.
Of note, the ability of GB models to yield better predictions than the PB model, in theory more accurate, is not a surprising result in light of previous works [67][68][69][70]. Since the GB OBC2 model accounted for the experimental relative affinities of the PPE and HNE complexes, this model was used in all further energy calculations, without an explicit mention of it hereinafter.

The effect of desolvation and protein-protein interactions on the formation of the studied complexes in solution
Polar-desolvation and screened electrostatic energies associated with the formation of the complexes in solution, in addition to the van der Waals and non-polar solvation energies, are shown in Table 2. These results show that ShPI-1 establishes stronger screened electrostatic interactions with HNE than ShPI-1/K13L (DDG sc el = 37.81 kcal/mol, Table 2). A similar behavior-although less pronounced-was obtained for the van der Waals and non-polar solvation energies (Table 2). However, the polar-desolvation penalty for the binding of ShPI-1 to HNE is larger than that of ShPI-1/K13L (ΔΔG ds = -42.63 kcal/mol, Table 2). This fact determines the more favorable ΔG eff value of the HNE:ShPI-1/K13L complex.
For the PPE complexes, a larger van der Waals energy contribution for the ShPI-1/K13L binding was predicted with respect to the 'in' conformation of the wild-type inhibitor (ΔΔE vw = -11.98 kcal/mol, Table 2). Conversely, the screened electrostatic interactions of the PPE:ShPI-1in complex are more favorable (DDG sc el = 7.88 kcal/mol, Table 2), but the global pair-wise interaction energy is still larger for the PPE:ShPI-1/K13L complex (DDG sc el þ DDE vw ¼ À4:10 kcal/mol). Note, however, that the relative contribution of pair-wise interactions does not account for the large ΔΔG calc value between these complexes, which mainly arises from their relative polar-desolvation penalty (ΔΔG ds = -15.59 kcal/mol, Table 2). On the other hand, there is a significant reduction of the polar-desolvation penalty in the PPE:ShP-1up complex with  Understanding the Affinity of PPE and HNE for ShPI-1 and ShPI-1/K13L respect to both the PPE:ShPI-1in an the PPE:ShPI-1/K13L complex (Table 2). Nonetheless, the van der Waals interactions and, more importantly, the screened electrostatic interactions are remarkably less favorable in the former complex (Table 2). Furthermore, the ΔG SA value is less favorable for the PPE:ShP-1up complex than for the other two PPE complexes (Table 2), indicating a decrease in the interface surface of the former. Overall, the weaker relative affinity of ShPI-1 for PPE is caused by the higher polar-desolvation penalty associated with the formation of the PPE:ShPI-1in complex, in addition to the less favorable pair-wise interactions at its interface when compared to those of the PPE:ShPI-1/K13L complex. An alternative 'up' conformation of the former complex is neither favorable due to the weak pair-wise interactions associated with its formation.
The impact of the K13L mutation on the polar-desolvation penalty. The formation of the PPE:ShPI-1in and the HNE:ShPI-1 complex occurs at expense of larger polar-desolvation penalties compared to those of both ShPI-1/K13L complexes ( Table 2). As observed, the greater affinity of the mutated variant for these enzymes is mostly determined by the ΔΔG ds values ( Table 2). To better understand the impact of the K13L mutation on the ΔΔG ds values, we performed a linear correlation analysis between the per-residue self-energy (DG rr GB ) and ΔSA res values for all positively-charged (Lys and Arg) and non-polar aliphatic (Leu, Ile, Val) residues in the studied complexes (Fig 9). According to our predictions, the molar energy required to desolvate a surface area of 1 Å 2 of positively-charged residues (8.20Á10 −2 kcalÁÅ -2 Ámol -1 ) is~3.60 times higher than that required for the non-polar aliphatic residues (2.28Á10 −2 kcalÁÅ -2 Ámol -1 ) (Fig 9). Thus, we inferred that the large DDG rr GB values (and, consequently, the ΔΔG ds values) associated with the K13L mutation are caused by the hydrophobicity differences between the Lys and Leu residues, together with the large desolvation underwent by both P1 residues upon their insertion into the S1 subsite, i.e., │ΔSA K13 │~230 Å 2 and │ΔSA L13 │~190 Å 2 , respectively (Fig 9).
Moreover, the results presented here indicate that the polar-desolvation penalty of K13 dramatically decreases in the 'up' conformation when compared to the 'in' conformation of the PPE:ShPI-1 complex, due to the reduction of the desolvation of this residue in the former conformation (Fig 9). Additionally, other residues at the PPE:ShPI-1up interface displayed lower desolvation penalties, especially those lying at the Pn and Sn sides of the inhibitor and the Linear Correlation between the Variation of Per-Residue Self-Energies and Solvent Accessible Surface Areas upon Complex Formation. The variation of per-residue self-energies (DG rr GB ) for the selected residues was obtained from the outputs of the pwEFED protocol, and correlated with corresponding opposite values of per-residue solvent-accessible surface area variation (-ΔSA res ) In turn, ΔSA res values were calculated from the ΔG SA(res) values present in the output of the prEFED protocol through the following equation ΔSA res = ΔG SA(res) /0.0072 (see Eq 1). The blue, orange and red dots on the left-hand panel correspond to K13 in the PPE:ShPI-1in, PPE:ShPI-1up and HNE:ShPI-1 complexes, respectively. On the right-hand panel, the red and blue dots correspond to L13 of ShPI-1/K13L in complex with HNE and PPE, respectively. The N-terminal residues I16 and V16 of HNE and PPE, respectively, were ruled out from the analysis, since the presence of a positively-charged amine group changes the overall hydrophobicity of these residues. The Pearson correlation coefficient for each linear fit is represented by r p . The linear-fit slopes represent an inverse measurement of the average hydrophobicity of each group of similar residues, i.e., the higher the slope value the lesser the residue hydrophobicity. Understanding the Affinity of PPE and HNE for ShPI-1 and ShPI-1/K13L enzyme, respectively (data not shown). Note, however, that the polar-desolvation decrease is concomitant to a reduction of the PPE:ShPI-1up pair-wise interaction energy (Table 2), which ultimately determines the low affinity of this complex.
Energy contributions of residues at the complex interfaces. The identification of warm-/ hot-spots at the interfaces of the studied complexes was first carried out using the prEFED protocol (S8, S9 and S10 Tables). This method predicted a relative free energy contribution of the P1 site residues in both HNE complexes significantly different from the ΔΔG calc value of the whole complexes (ΔΔG sc(L13-K13) = -7.52 kcal/mol and ΔΔG calc = -0.79 kcal/mol, see S8 Table  and Table 1A, respectively). Thus, we performed the CAS protocol for all the residues previously identified as warm/hot-spots by the prEFED approach. Of note, we excluded from this procedure the Cys, Pro, Gly and Ala residues, for which the mutation by Ala is either structurally forbidden or meaningless.
The predictions of CAS and prEFED were linearly correlated (r p = 0.90) and similarly ranked the residues according to their energy contributions (r s = 0.87) (Fig 10A), in accordance with previous reports [60]. As observed, the points corresponding to D226 and K13 of the HNE:ShPI-1 complex are outliers of the linear fit, since their elimination from the statistical analysis improved the values of the correlation coefficients ( Fig 10B). This fact proved that the energy predictions of CAS and prEFED significantly differed for both residues. Moreover, only CAS yielded a relative free energy contribution for K13 and L13 at the P1 site of both inhibitors consistent with the ΔΔG calc value (ΔΔG (L13-K13) = -0.46 kcal/mol). This discrepancy is likely to be caused by the well-known proneness of GB models to errors in the calculation of individual DG ij GB values, which nearly cancel during the calculation of total values ( X i;j DG ij GB ) [71][72][73]. Since CAS predictions are based on the ΔG eff values for the native and the mutated complex, it becomes apparent that they are less influenced by internal errors in the DG ij GB calculations than those obtained through prEFED. Overall, we concluded that CAS was more robust than prE-FED, especially when predicting the energy contributions of residues involved in salt-bridge interactions, like D226 and K13 of the HNE:ShPI-1 complex.
From the CAS results, three residues of the S1 subsite, i.e., H57, Q/F192 and D194, were identified as the most important common/equivalent hot-spots of PPE and HNE (ΔΔG -4.0 kcal/ mol) (Fig 11). Other S1 residues conserved in both elastases, i.e., S195, S214, F215 and V216, were also predicted as warm-/hot-spots (Fig 11), thereby reinforcing the importance of this subsite in the interaction with the inhibitors. Additionally, the residue D226, also located within the S1 subsite of HNE, had the largest energy contribution to the formation of the complex with the Understanding the Affinity of PPE and HNE for ShPI-1 and ShPI-1/ K13L  Fig 11. Per-Residue Energy Contributions to the Formation of the Studied Complexes. Residues identified as warm/hot-spots by prEFED were selected for CAS, except for Cys, Pro, Ala and Gly. The vertical dashed lines separate the enzyme's residues from those of the inhibitor. X13 stands for L13 or K13 depending on the complex. ΔΔG = ΔG eff (native complex)-ΔG eff (mutated complex), therefore, a negative ΔΔG value indicates a favorable energy contribution of the mutated residue to the complex formation. Conversely, a positive ΔΔG value indicates and unfavorable contribution of the mutated residue. Error bars were not included since standard errors were always lesser than 5% of the mean values. A structural representation of the warm-/hot-spot residues at the interface of each complex is shown. Residues were colored according to their energy contribution to the complex formation, see colorgradient bars. doi:10.1371/journal.pone.0137787.g011 Understanding the Affinity of PPE and HNE for ShPI-1 and ShPI-1/K13L wild-type inhibitor (Fig 11). In fact, the mutated enzyme ENH/D226A would be unable to bind ShPI-1 due to an increase of~16 kcal/mol in the ΔG eff value (Fig 11). In contrast, the equivalent residue of PPE, T226, is not involved in the hypothetical interaction of this enzyme with ShPI-1. On the other hand, the residue at position 226 of both elastases is irrelevant for their interaction with ShPI-1/K13L (Fig 11). Hence, these results suggest that the mutated enzymes PPE/T226A and HNE/D226A display a PPE-like behavior against both inhibitor variants. For ShPI-1 and ShPI-1/K13L, highly-favorable energy contributions were predicted for the P1 site residues in all complexes except for both conformations of the PPE:ShPI-1 complex (Fig 11), which is not actually formed [7]. This result is consistent with the crucial importance of the P1 site residue in the interaction of canonical inhibitors with target proteases [13]. Interestingly, the energy contribution of K13 is less unfavorable in the 'up' conformation than in the 'in' conformation, in accordance with its lower polar-desolvation penalty when adopting the former conformation. Finally, warm-/hot-spots outside the S1:P1 interfaces of the complexes were also identified (see S2 Text). Remarkably, residue E44(P31') lying outside the binding loops of the inhibitors was identified as a hot-spot at the interfaces of the HNE complexes (Fig 11).
To clarify the origin of the energetic contribution of these warm/hot-spots residues in terms of their pair-wise interactions with neighboring residues, the pwEFED protocol was employed. Residues D226 of HNE and K13 of ShPI-1 were predicted as the pair of interacting residues with the largest energy contribution to the HNE:ShPI-1 complex formation (ΔG D226, K13 = -26.87 kcal/mol). This result confirms the importance of the hydrogen bonds and salt bridge formed between the previous residues ( Fig 4B) for the binding of the wild-type inhibitor to HNE. In agreement with the results obtained through the prEFED and CAS approaches, the negligible energy contribution of D226 of HNE to the binding of ShPI-1/K13L arises from its weak pair-wise interaction with L13 at the P1 site (ΔG D226, K13 = -0.37 kcal/mol). A similar result explains the slight contribution of T226 of PPE in the interaction of this enzyme with both inhibitor variants (ΔG T226, K13in = -0.07 kcal/mol, ΔG T226, K13up = 0, ΔG T226, L13 = -0.61 kcal/mol). Interestingly, we also observed a favorable pair-wise interaction between E44 (P31') of either ShPI-1 or ShPI-1/K13L and R36 of HNE (ΔG R36, E44 = -3.65 kcal/mol and -5.30 kcal/ mol, respectively) due to the hydrogen bond and salt bridge formed between these residues (S3 and S4 Tables).
Interestingly, we also predicted that the total energy contribution of L13 due to pair-wise interactions with the PPE residues (-33.02 kcal/mol) is less than that of K13 in the PPE:ShPI-1in complex (-36.00 kcal/mol). This counterintuitive result highlights the impossibility of correlating the strength of pair-wise interactions with the complex affinity without taking into account the polar-desolvation penalty effects. Conversely, the pair-wise energy contribution of K13 in the PPE:ShPI-1up complex (-25.72 kcal/mol) is smaller than those of L13 and K13 in the PPE:ShPI-1/K13L and PPE:ShPI-1in complexes, respectively, in accordance with the reduced number of contacts established by K13 in the 'up' conformation ( Fig 6B). It is noteworthy that the K13L mutation at the P1 site caused changes in the pair-wise interaction energies between residues outside the S1:P1 interface (data not shown). The occurrence of such changes spread over the entire complex interfaces, but mostly at the P2-P7 sites of the inhibitor and the corresponding subsites of the enzyme, helps understand why the PPE:ShPI-1/K13L complex possesses a larger pair-wise energy contribution compared to the PPE:ShPI-1in complex, even though L13 establishes weaker pair-wise interactions than K13. For the HNE complexes, pair-wise interaction energies of -33.00 kcal/mol and -68.98 kcal/mol were predicted for L13 and K13, respectively. However, the energy contribution of the latter residue is greatly diminished by its large polar-desolvation penalty (Fig 9). In addition, changes in the pair-wise interaction energies outside the S1:P1 interface, especially comprising the Sn' and Pn' sides of HNE and both inhibitor variants, respectively, (data not shown) caused the more favorable van der Waals energy of the HNE:ShPI-1/K13L complex with respect to the HNE:ShPI-1 complex (see Table 2).

Discussion
Position 226 determines the differential specificity of HNE and PPE for ShPI-1 The structural and energetic analyses performed in this work pointed out the critical role of position 226 of HNE and PPE in determining their ability to interact or not with the wild-type inhibitor ShPI-1. According to 3D structure model of the HNE:ShPI-1 complex, the residue D226 at the S1 subsite of HNE forms relatively stable hydrogen bonds and a salt bridge with the basic residue K13 at the P1 site of ShPI-1 (Fig 4B). This interaction is essential to overcome the high polar-desolvation penalty of the K13 residue (Fig 9). Furthermore, these predictions suggested that the enzyme variant ENH/D226A would be unable to bind ShPI-1 (Fig 11), thereby proposing an experimental approach to demonstrate the crucial role of the D226-K13 interaction in the complex formation.
The presence of D226 at the bottom of the S1 subsite of HNE would also explain the interaction of this enzyme with other canonical inhibitors bearing Lys or Arg at the P1 site, such as BPTI and the Kazal-type inhibitor CmpI-II, respectively [20,74]. Indeed, the occurrence of close contacts ( 4 Å) between D226 and R12(P1) of the latter has been proposed before based on the analysis of a 3D-structure model of the HNE:CmPI-II complex [74]. Furthermore, in accordance with our predictions based on a GB model, a previous study also assessed the high desolvation penalties of ligands with positively-charged moieties at the P1 site upon their binding to thrombin, using quantum chemical calculations of solvation energies. A similar role to that of D226 was also suggested for residue D189 of thrombin in the compensation of the desolvation penalties of positively-charged moieties at the P1 site [75].
Unlike HNE, PPE has a Thr residue at position 226 that is unable to interact with K13 of the wild-type ShPI-1 (Fig 6). The previous fact may preclude the formation of the PPE:ShPI-1in complex, since the high polar-desolvation penalty associated with the insertion of K13 into the S1 subsite of PPE is not compensated by strong pair-wise interactions. An alternative 'up' conformation of K13 was sampled during the 125 ns MD simulation of the PPE:ShPI-1 complex (Fig 5). Interestingly, a similar conformation has been observed in the crystal structures of ShPI-1 and BPTI in complexes with bovine chymotrypsin (PDB: 3T62 and 1CBW, respectively), another SP with a hydrophobic S1 subsite [7,15]. Of note, S217 at the entrance of the S1 subsite of chymotrypsin helps stabilize the 'up' conformation of a Lys residue at the P1 site through the formation of a hydrogen bond [7,15]. In contrast, our predictions revealed the lack of strong interactions of K13 of ShPI-1 with the residues at the entrance of the S1 subsite of PPE. In particular, the absence of a hydrogen bond between K13(NZ) and S217(O) in the PPE:ShPI-1up complex (S6 Table), was corroborated, which is consistent with previous suggestions [7].
On the other hand, the results obtained from the MD simulations of the disrupted complexes also revealed the importance of the long-ranged electrostatic interaction between D226 of HNE and K13 of ShPI-1 in the binding process. In fact, upon an initial interface disruption in the HNE:ShPI-1 complex, the distance between both residues decreased from~5.7 Å to an equilibrium distance of~2.7 Å during the simulation time (Fig 7), which suggests that the D226-K13 attraction may be an essential driven force for the insertion of K13 side-chain into the enzyme S1 subsite. Conversely, the absence of long-range interactions between K13 and T226 of PPE is likely to preclude the formation of the PPE:ShPI-1in complex (Fig 7).
Our results strongly suggest the importance of the residue at position 226 of HNE and PPE in determining whether a positively-charged group of the inhibitor can be accommodated into the S1 subsite of these enzymes. In addition, the structural characteristics of the S1 subsite of PPE hamper the stabilization of inhibitor's Lys residues at the entrance of this subsite. We suggest that these facts are the underlying causes of the experimentally demonstrated large difference in the affinity of ShPI-1 for HNE and PPE.
The lower polar-desolvation penalty contributes to the higher affinity of ShPI-1/K13L for both elastasesd The higher affinity of ShPI-1/K13L for PPE and HNE compared to that of ShPI-1 depends to a great extent on the lower polar-desolvation penalties associated with the binding of the former inhibitor to both enzymes ( Table 2). In fact, this is the only energy component explaining the higher affinity of ShPI-1/K13L for HNE ( Table 2). On the other hand, the PPE:ShPI-1/K13L complex has both stronger pair-wise interactions and a less unfavorable polar-desolvation penalty than the PPE:ShPI-1in complex. Interestingly, the stronger pair-wise van der Waals interactions of ShPI-1/K13L with PPE compared to ShPI-1 are not directly caused by its more hydrophobic P1 site residue, but rather by other neighboring residues at the Pn side of its primary binding loop, e.g., V9(P5), G10(P4) and, especially, R11(P3).
The well-known preference of elastases for aliphatic residues at the P1 site [23,76] has been usually attributed to the relative hydrophobicity of their S1 subsites [22] (S1 Fig). Our results showed that this hydrophobicity should be understood as the lack of residues at key positions capable of establishing pair-wise interactions with polar P1 site residues, strong enough to overcome the highly unfavorable polar-desolvation penalties of the latter. Accordingly, the hydrophobicity of a binding site cannot be viewed as an average property mainly determined by the nature of its most abundant residues. In fact, the sole difference at position 226 has a great impact on the binding of HNE and PPE to PIs bearing basic residues at the P1 site; although, the S1 subsites of both enzymes have been commonly described as hydrophobic [22]. P3, P2' and P3' sites as important positions for modulating the inhibitor specificity and affinity for elastases As expected for a canonical serine protease inhibitor, the warm-/hot-spots residues of ShPI-1 and ShPI-1/K13L in their interaction with both elastases were identified mainly within the primary and secondary binding loops [13]. Interestingly, from all the previously identified residues, those at the P3, P1, P2' and P3'sites were predicted to interact with regions of HNE and PPE displaying high variability in their amino acid composition (S2 Table). Considering such differences, amino acid substitution at the previous sites could modify the specificity of the inhibitor toward HNE or PPE. This has been experimentally confirmed at least for the P1 site of ShPI-1 and suggested for the P3 site of this inhibitor based on a previous structural analysis [7]. Moreover, our predictions show that R11 at P3 of ShPI-1 has more favorable energy contribution in the PPE complexes, mainly due to the presence of the negatively-charged residue D98 within the S3 subsite of PPE (N99A in HNE) (see S1 Fig and Fig 11). Thus, mutations at this position may enhance the specificity of ShPI-1/K13L for HNE.
Besides the well-known involvement of the P1 site residue in the binding to all SPs, including elastases, residues Y15(P2') and F16(P3') were predicted to have large energy contributions, comparable to those of the P1 site in the three existing complexes analyzed here (Fig 11). This result underscores the critical role of the interactions at the S2':P2' and S3':P3' interfaces in elastase inhibition, although a subordinate role of the interactions at the Sn' and Pn' sides with respect to those at the Sn and Pn sides of SPs and their substrate-like inhibitors, respectively, has been previously referred [77]. It is also likely that both elastases display a preference for aromatic residues at P2' and P3' sites. Hence, unlike the S3 subsite, the different residue composition between the respective S2' and S3' subsites is not expected to promote differences in the specificities of HNE and PPE for inhibitor residues at P2' and P3' sites. The preference for aromatic residues at P2' and P3' sites could explain the higher affinity of ShPI-1 for HNE (K i = 1.3Á10 −9 M) compared to BPTI (K i = 3.5Á10 −6 M) [7,20], since the latter possesses non-aromatic Arg and Ile residues at these subsites.
E44: A hot-spot located outside the inhibitor's standard binding loops Apart from the residues of the primary and secondary binding loops, our predictions suggest a hot-spot residue outside both loops, i.e., E44(P31'), which is involved in a hydrogen bond and a salt bridge with R36 of HNE in the complexes with ShPI-1 and ShPI-1/K13L (S3 and S4 Tables). The polar interactions between these residues constitute the structural basis of their energy contribution to the complex formation (Fig 11).
To our knowledge, this is the first work reporting the putative involvement of a residue so distal from the binding loops of a BPTI-Kunitz domain in the interaction with a target protease. It is noteworthy that a previous report has suggested the impact of residues near E44 (K46 in BPTI), since mutations within the segment 39-RAKR-42, contiguous to the secondary binding loop of BPTI, influence the inhibitory activity of this inhibitor against HNE [78]. However, the influence of K46(P31') on the interaction of BPTI with HNE has not been explored to date. On the other hand, according to our predictions, E44 has an unfavorable contribution to the formation of the complex with PPE (Fig 11). Hence, we hypothesize that amino acid substitution at this position would contribute to the design of ShPI-1/K13L variants with increased specificity for PPE.

Conclusions
In this work the structural and energetic analyses of the interaction of ShPI-1 and ShPI-1/K13L with PPE and HNE were performed. Residue D226 was identified as structural hallmark of HNE for the interaction with inhibitors with basic chains at the P1 site. The lower polar-desolvation penalties associated with the binding of the ShPI-1/K13L to both elastases are also suggested to constitute the main energy component contributing to its larger affinity for these enzymes compared to that of the wild-type inhibitor ShPI-1. In addition, we carried out an extensive study of the interactions and the energy contributions of the residues at the complex interfaces. This analysis confirmed the importance of residues of the primary and secondary binding loops of BPTI-Kunitz inhibitors in their interaction with SPs and, especially, that of the P1 site residue. In addition to this site, residues at the P3, P2' and P3' were predicted as important positions to modulate the specificity and/or affinity of the inhibitors for elastases. We also suggested for the first time that E44, a residue located at a distal position from the binding loops, has a significant energy contribution to the interaction with HNE. Overall, the methodology and results of this work will be valuable for the design of new ShPI-1 variants more potent and/or specific for either HNE or PPE.
Supporting Information S1 Fig. Structural Alignment of PPE and HNE Shown at the Level of Their Respective Primary Sequences. The sequence identity between both proteins is~39%. Asterisks have been placed under the conserved positions and black dots, under the residues belonging to the S1 subsite of each elastase. Chymotrypsinogen residue numbering has been adopted [22]. Residues depicted in red belong to the catalytic triad of SPs. Additionally; yellow rectangles have been used to highlight the conserved Cys residues. The structural alignment was carried out with Modeller v9.5. . RL, R and L stand for the complex, the receptor and the ligand, respectively. The blue box represents the solvent. ΔE el is the electrostatic energy variation upon complex formation in vacuum. G GB(X) represents the polar-solvation free energy of molecule X, where X stands for R, L or RL. Similarly, G ii GBðXÞ represents the polar-solvation free energy of atom i belonging to molecule X. The partial charges of all atoms of R and L except for atom i (orange dot) were set to zero. The inner space of R and L is shown in white to suggest the absence of intra-solute electrostatic interactions. The iteration of the thermodynamic cycle for every particle i of the solute molecules leads to the calculation of ΔG ds as indicated in the figure. Note, however, that two terms standing for cross-energies (DG ij GB ) between particles i and j within the same solute molecule (R or L) must be added to the self-energies (DG ii GB ). (TIF)  Table. Comparison of the van der Waals Contacts and Hydrogen Bonds at the Interfaces of the PPE:ShPI-1/K13L Crystal and Representative Structures. Van der Waals contacts were determined with a cutoff radius of 4 Å. For hydrogen bonds, the geometric constraints were i) a distance 3.5 Å between the donor and the acceptor and ii) an acceptor-donorhydrogen angle 30°. The names of the donor and acceptor atoms, as well as the donor-acceptor distance and the hydrogen bond occupancy during the productive MD simulation are shown in red. The names of PPE residues involved in Van der Waals contacts conserved in both complexes are shown in bold style. (DOCX) S2  Table. Main Per-Residue Energy Contributions Calculated by prEFED for the HNE:ShPI-1 and HNE:ShPI-1/K13L Complexes. Only those residues for which |ΔG res |!1.0 kcal/mol and/or | ΔG sc |!0.5 kcal/mol at least in one of the complexes are shown here. Note that E44 was included because of its interaction with R36, even though it does not fulfill the previous conditions. (DOCX) S9 Table. Main Per-Residue Energy Contributions Calculated by prEFED for the PPE: ShPI-1in and PPE:ShPI-1/K13L Complexes. Only those residues for which |ΔG res |!1.0 kcal/ mol and/or |ΔG sc |!0.5 kcal/mol at least in one of the complexes are shown here. Note that T226 was included to compare its energy contribution wih that of D226 of HNE, even though it does not fulfill the previous conditions. (DOCX) S10 Table. Main Per-Residue Energy Contributions Calculated by prEFED for the PPE: ShPI-1up and PPE:ShPI-1/K13L Complexes. Only those residues for which |ΔG res |!1.0 kcal/ mol and/or |ΔG sc |!0.5 kcal/mol at least in one of the complexes are shown here. Note that T226 was included to compare its energy contribution wih that of D226 of HNE, even though it does not fulfill the previous conditions. (DOCX) S1 Text. Analysis of the Polar Interactions Outside the S1:P1 Interfaces in the Three Existing Complexes. (DOCX) S2 Text. Per-Residue Energy Contributions outside the S1:P1 Interfaces of the Studies Complexes. (DOCX)