How does heparin prevent the pH inactivation of cathepsin B? Allosteric mechanism elucidated by docking and molecular dynamics
© Costa et al. 2010
Published: 22 December 2010
Skip to main content
© Costa et al. 2010
Published: 22 December 2010
Cathepsin B (catB) is a promising target for anti-cancer drug design due to its implication in several steps of tumorigenesis. catB activity and inhibition are pH-dependent, making it difficult to identify efficient inhibitor candidates for clinical trials. In addition it is known that heparin binding stabilizes the enzyme in alkaline conditions. However, the molecular mechanism of stabilization is not well understood, indicating the need for more detailed structural and dynamic studies in order to clarify the influence of pH and heparin binding on catB stability.
Our pKa calculations of catB titratable residues revealed distinct protonation states under different pH conditions for six key residues, of which four lie in the crucial interdomain interface. This implies changes in the overall charge distribution at the catB surface, as revealed by calculation of the electrostatic potential. We identified two basic surface regions as possible heparin binding sites, which were confirmed by docking calculations. Molecular dynamics (MD) of both apo catB and catB-heparin complexes were performed using protonation states for catB residues corresponding to the relevant acidic or alkaline conditions. The MD of apo catB at pH 5.5 was very stable, and presented the highest number and occupancy of hydrogen bonds within the inter-domain interface. In contrast, under alkaline conditions the enzyme's overall flexibility was increased: interactions between active site residues were lost, helical content decreased, and domain separation was observed as well as high-amplitude motions of the occluding loop – a main target of drug design studies. Essential dynamics analysis revealed that heparin binding modulates large amplitude motions promoting rearrangement of contacts between catB domains, thus favoring the maintenance of helical content as well as active site stability.
The results of our study contribute to unraveling the molecular events involved in catB inactivation in alkaline pH, highlighting the fact that protonation changes of few residues can alter the overall dynamics of an enzyme. Moreover, we propose an allosteric role for heparin in the regulation of catB stability in such a manner that the restriction of enzyme flexibility would allow the establishment of stronger contacts and thus the maintenance of overall structure.
Cathepsin B (EC 22.214.171.124) (catB) is one of the most well-characterized cysteine proteases, belonging to the clan CA (papain superfamily). In humans, its physiological role is implicated in bone resorption, antigen processing and protein turnover . However, catB also participates in pathological processes such as cardiovascular disturbances , parasitic infections , Alzheimer's disease , osteoarthritis , tumor invasion and metastasis development [6, 7]. Its main roles in cancer are i) its activity in directly cleaving extracellular matrix (ECM) components, ii) its activation of other ECM degrading proteases, which promotes tumor cell invasion into the surrounding tissue and bloodstream escape , and iii) stimulating angiogenesis which provides increased nutrients and oxygen supplies to cancer cells . Thus, catB regulates several crucial steps in tumorigenesis and is a promising target for anti-cancer drug design .
Currently the most potent and selective structure-based designed compounds available are derived from E-64 targeting the unusual occluding loop present in the catB 3D-structure . However, enzymatic assays have shown that these inhibitors are strongly pH dependent as their optimal binding affinities are considerably diminished under neutral/alkaline conditions . Since catB can be found in several cellular compartments with distinct pH values, these inhibitors are not effective in vivo. When catB is within the lysosomal or endosomal vesicles it confronts acidic conditions, in contrast to the neutral/alkaline environment when it is attached to membranes (mainly in caveolar microdomains) [17, 18] or secreted in the ECM . Although catB is rapidly inactivated under alkaline pH conditions, it was reported that membrane-associated forms are resistant to this process . This peculiarity is believed to occur due to its interaction with heparan sulfate glycosaminoglycans (GAGs) on the cell surface . This polysaccharide, structurally related to heparin, acts on the ECM as a co-receptor for several molecules such as growth factors and proteases . Interaction between catB and heparin-like GAGs was already shown to prevent the enzyme's inactivation in high pH . The main reported outcome was the maintenance of catB helical content in the presence of heparin at high pH, which was observed for papain as well [22, 23]. Nevertheless, it has not been possible to precisely define the catB/heparin interaction sites and the molecular mechanism responsible for this protective effect.
Structural and molecular modelling studies can give insight into the molecular events concerning the modification of catB activity by pH changes and heparin binding. Some attempts have already been made towards designing new catB inhibitors using molecular modeling techniques [24, 25], taking into account dynamical aspects of binding. However, these studies did not evaluate the influence of pH nor that of heparin binding on the modulation of the enzyme intrinsic flexibility.
In order to understand the heparin protective effect over catB structure we performed molecular dynamics (MD) simulations using an approach in which different pH conditions were taken into account by considering different protonation states of titratable groups on the protein surface. Further, docking analyses resulted in the identification of two heparin binding sites on catB structure. The MD calculations confirmed an increase of the overall catB flexibility and the loss of stability of the apo catB inter-domain interface in alkaline pH. The destabilization and the increased flexibility, notably in the occluding loop, were prevented by interaction with heparin, again in agreement with experimental data. We observed a role of active site residues in enzyme stabilization and in maintaining the helical content, and we propose an allosteric mechanism for the stabilizing effect promoted by GAG interaction. Taken together, our data provides an improved understanding of the molecular mechanisms responsible for both pH-induced inactivation and protection against inactivation by heparin binding.
Prediction of pKa's for protein ionizable residues is an important tool for understanding features and catalytic mechanisms of pH-dependent enzymes . We applied the PROPKA program to estimate pKa values in the catB structure and to determine the most probable corresponding protonation states of the enzyme under acidic and alkaline conditions. Two different conditions were studied: acidic (pH 5.5), and alkaline (pH 8.0), which allow comparison to available experimental data .
It is known that catalytic residues often present unusual pKa values compared to those of free amino acids in solution . Accordingly, our results for catB predicted a pKa of 2.14 for the catalytic C29, in contrast to 8.0 expected in solution. Previous work on the papain catalytic cysteine (structurally related to C29 in catB) showed pKa values around 2.5 to 3.5 , in accordance to our estimation for C29.
Differentially protonated residues in catB and pKa values
Average values during 40ns MD
3.7 ± 1.5
5.4 ± 0.5
3.4 ± 1.2
7.1 ± 0.2
7.3 ± 0.4
6.6 ± 0.8
8.1 ± 0.7
7.4 ± 0.2
6.5 ± 0.4
7.1 ± 0.2
7.4 ± 0.7
8.1 ± 2
4.5 ± 1.7
As expected, most of the residues with differing protonation states were histidines, with a pKa in solution of 6.5 and an ionization state that is very susceptible to pH changes in the physiological range. H110 is a key residue for stabilization of the occluding loop and is also crucial for the exo-proteolytic activity of catB, since along with H111 it anchors the carboxyl terminus of substrates. Interestingly, we observed an pKa of 7.79 for the catalytic H199, which implies that this residue is protonated only in acidic conditions. This result is correspondent to the measured pKa of H159 in papain . In papain, it was previously reported that the ionization of the catalytic C25 and H159 residues are coupled . There the deprotonation of H159 shifts the pKa of C25 from ~3.3 to 7.6, while the neutralization of C25 decreases the pKa of H159 from ~8.5 to 4.3. Due to the direct relationship between catB and papain we assigned the protonation states of these residues according to this proposed mechanism.
The pKa values predicted with PROPKA thus appear to be sensitive to local microenvironment changes around the ionizable residues. In order to verify if the applied protonation profiles accurately represent the distinct pH conditions simulated, we collected 1000 snapshots (one at each 400 ps) from each molecular dynamics (MD) trajectory in order to perform pKa calculations. Table 1 shows the average pKa values during the MD simulations. From both pH conditions we obtained values consistent with the protonation states predicted from the crystal structure, thus confirming their validity. In this Table we have also highlighted the pKa shift in the catalytic residues C29 and H199 as observed for papain.
In MD simulations of the bound enzyme at pH 5.5, the RMS fluctuation profile was seen to be similar to that of the apo form (Fig. 4B). Hence, catB is stable under acidic conditions independent of heparin binding. In contrast, under alkaline conditions heparin-bound catB exhibited much smaller backbone fluctuations (2 Å) than the very flexible apo form in these conditions. Moreover, we observed similar results when simulations were started using different initial conformations and also with heparin docked at the other binding site (L domain) (Fig. 4B). Functional regions including the occluding loop and the active site were also stabilized, showing that heparin binding at the predicted sites promotes a global stabilizing effect.
In Fig. 4C, we show the time-evolution of the number of distinct clusters of catB backbone conformations throughout the different MD simulations. Here each cluster contains conformations within an RMSD of 1 Å of the cluster center. If during a simulation the system visits numerous clusters, the RMSD from the starting structure will be very high, indicating conformational instability in the simulation. In contrast, if a simulation visits only a few (correspondingly larger) clusters, it can be considered as being more conformationally stable. Indeed, we found few clusters during the simulations of apo catB and the catB-heparin complex (4 and 2 clusters, respectively) under protonation conditions corresponding to pH 5.5 (Fig. 4C). Furthermore, in acidic conditions the backbone RMSD distributions, similar to the RMS fluctuation profiles, did not change significantly on binding heparin, presenting deviations in both cases centered at 2.2 Å. However, in the heparin-bound simulation in acidic conditions, the RMSD distribution began to shift towards a conformational state far from the native structure. (Fig. 4D).
Considering that catB loses stability and activity in alkaline conditions , in the systems simulated under protonation conditions corresponding to high pH we expected to observe a higher number of backbone clusters. This expectation was confirmed: we observed 66 clusters in the apo catB simulation (Fig. 4C). In addition, under these conditions the catB structure does not stabilize during the simulation, with a continuous increase of backbone RMSD values throughout the 40 ns production period. A similar profile was observed in the replica simulation and also when starting with another conformation (a separate PDB structure, see Methods), in which we found 47 clusters. The RMS fluctuations of both replica simulations, plus the simulation performed with the alternative catB structure (Additional file 1), showed similar results to the apo enzyme in alkaline conditions (Fig. 4A).
The RMSD distribution plot shows a clear shift dividing into two populations: one around 2 Å and the other at higher RMSD values. (Fig. 4D). On the other hand, in the simulation with heparin under alkaline conditions we observed fewer clusters comparable to acidic conditions (Fig. 4C). This stabilization afforded by heparin binding is also reflected in the RMSD distribution presenting a narrowed profile centered at small RMSD values (around 1.7 Å) (Fig. 4D).
The results shown above suggest that catB as a whole and in particular the critical domain interface and occluding loop are stabilized by heparin binding to a relatively distant site in the R domain. From basic thermodynamic principles, preferential binding of ligand to the native state will generally confer stability to a macromolecule in what has been described as an allosteric mechanism, without necessity for direct interaction between the ligand binding site and the active site residues, and has been extensively described in phenomenological terms . The structural information available for catB, coupled with macromolecular simulations, provides a means for investigating details of such a mechanism directly.
We applied an essential dynamics analysis to examine the most relevant global macromolecular motions occurring during the simulations. Such analyses have been largely applied in the understanding of biological functions in proteins since they provide an evaluation of large-scale movements that are often related to domain motions and conformational changes . We thus diagonalized the atomic positional covariance matrix to obtain the eigenvectors and corresponding eigenvalues. Selecting the first five (largest amplitude) principal components, we could recover around 60% of the total motions for apo catB and around 56% for heparin-bound catB. Our analysis focused on the systems simulated under alkaline conditions in order to address how heparin influences catB flexibility and prevents domain separation as shown above. We point out that certain analyses involving averaged quantities such as the essential dynamics and clustering approaches can be best interpreted for stable systems; the precise results for the apo catB under alkaline conditions will to some extent depend on the starting conformation. However, even under these conditions the overall results obtained were similar for replicated runs.
By inspecting other motions corresponding to the first component (PC1) in the R-domain complexes, we identified local changes in the binding region. The residues close to the polysaccharide (L1, A3 K156 and R233) adopted a new conformation. To verify that these motions would not lead to loss of binding interaction between catB and heparin, we measured the distance between the center of mass of heparin and the positive site during MD (Additional file 2) and found stable behavior: heparin remained bound during the entire simulation. Remarkably, however, this local motion is coupled to movements far from the heparin binding site. We observed that the residues 58 to 75 in the L-domain moved towards the interdomain interface (Fig. 6C). This motion certainly contributes to the stabilization of the overall structure, thus helping explain the results described above (Fig 4 and Fig 5). We did not find equivalent behavior in the L-domain complexes, which led us to investigate more deeply the correlated motions of catB in the apo form and in complex with heparin.
The analysis of cross correlation coefficients between pairs of residues is a useful method to identify correlated collective motions [37, 38]. The correlation matrix represents the linear correlation between pairs of C-alpha atoms as they move about their average positions during dynamics. Positive correlations are related to motions occurring in the same direction whereas negative correlations indicate motion in opposite directions. We compared the correlation matrix of apo catB and of the complexes. In the apo form we observed positive intra-domain correlations but negative inter-domain correlations (Additional file 3). This pattern is not observed in the complexes, in which we identified a more diffuse pattern of correlations regardless of the binding site occupied. This result shows that heparin acts mainly as a stabilizing element, and its binding seems to restrict the anti-correlated collective motions that would be responsible for domain separation.
Concerning the occluding loop, we observed that the opening-closing movement is clearly represented by the two first components in apo catB. On the other hand, upon heparin binding these motions were not observed, independent of the binding site occupied. In addition, the first PC reveals the increasing distance between C29 and the other catalytic residues in apo catB, while in the complex this motion is absent. Therefore, heparin stabilizes catB motions in its functional regions. These results are consistent with experimental findings, which, although they revealed the implications of binding on the activity of the enzyme, did not provide a structural view of the phenomenon . Heparin binding, despite the distance separating it from active site residues, prevents their disorganization as would otherwise be produced under alkaline conditions. The principal components obtained from the essential dynamics analysis of our MD simulations suggest that global concerted movements in the macromolecule leading to domain separation are suppressed by heparin binding, thus helping explain an allosteric role for heparin in this biological system.
Loss of helical content was reported as the main effect related to alkaline-pH-induced inactivation of catB . In the same study, it was reported that heparin prevented the loss of helical structure, which led to the postulate that the maintenance of activity in the catB-heparin complexes is strictly related to this structural aspect.
It is discussed in the literature that papain-like enzymes, which possess a similar domain organization, present structural characteristics that confer some rigidity to the active site region . The most important element supporting this characteristic is the network of interactions in the interdomain interface. We verified the stability of the interaction between catalytic residues C29-H199 and found that disorganization of the main alpha helix in apo catB under alkaline conditions is coupled with the separation of the catalytic residues (Additional file 5). These events are strongly correlated with the overall domain motions observed. The binding of heparin at a substantial distance form the interface nevertheless stabilizes the interdomain contacts maintaining the structuring of the active site, and would explain the results obtained in biochemical assays .
The occluding loop controls access to the active site and also confers the exo-proteolytic activity of the enzyme. In our simulations, we verified the open-close mechanism of this region and found that in acidic conditions a stable closed state is maintained regardless of heparin binding. However, at pH 8.0 the occluding loop exhibited high flexibility in the absence of heparin (Fig. 7B) as expected after RMSF analysis (Fig 4A). Nevertheless, under alkaline conditions, when heparin was bound the occluding loop exhibited a more stable conformation (Fig. 7B). In this conformation, besides closing the active site the occluding loop interacts with the R-domain in a distinct fashion from that observed in apo catB (Additional file 4). The occurrence of this conformation was also independent of the positive site occupied by heparin since the RMSD between the different complexes is around 1.4 Å. On the other hand, the loop conformations observed in the apo form simulations exhibit higher RMSD values (around 4.5 Å) compared to the complexes.
This result allows a discussion of experimental findings in which it was shown that heparin inhibits the exo-proteolytic activity of catB . In that study, the authors proposed a competitive mechanism of inhibition due to the possible interaction of heparin with H110 and H111, which are the most important residues involved in accommodation of the negatively-charged C-terminal substrate carboxylate groups. From our analysis we conclude that the mechanism of inhibition promoted by heparin appears to be related also to the induction of conformational changes in catB, mainly in the occluding loop region, which adopts a distinct stable conformation that may impair proper binding of small substrates.
However, we cannot exclude the possibility that large heparin polymers may also bind to the occluding loop, thus inhibiting the exo-proteolytic activity of the enzyme, especially if the concentration of heparin increases as shown in ref. .
Occupancy of hydrogen bonds between interacting interface residues
Percentage of occupancy of hydrogen bonds between interface residues
apo catB (pH 5,5)
Apo catB (pH 8,0)
catB + hep (pH 5,5)
catB + hep (pH 8,0)
W30 – E171
W80 – E171
D81 – S150
R41 – P169
Q23 – S220
These results are correlated with the distinct loop conformation observed in the heparin bound complexes (Additional file 4). In addition, as shown, heparin binding at the R-domain induces a conformational change that leads to a more compact structure (Fig 6C). Further, heparin restricts collective motions involved in domain separation, which allows the stabilization of a distinct pattern of interactions formed in the interface under alkaline conditions.
Heparin-protein interactions are known to regulate several biological processes including protease regulation , growth factor activity , macromolecular assembly  and viral mechanisms . In this paper we unraveled the molecular mechanism of heparin protection against pH-induced inactivation of catB. This phenomenon was previously demonstrated through biochemical assays without a full understanding of the process because of the lack of structural information concerning heparin binding. Herein we were able to mimic the different pH conditions by performing MD simulations of catB with different protonation states according to pKa calculations. We confirmed experimental findings showing that under acidic conditions catB is stable. In contrast, at alkaline pH six residues are deprotonated, increasing the negative charge of the interdomain interface. This leads to charge repulsion and subsequent separation of the L and R domain. In addition, we observed the effect of the alkaline conditions in terms of destabilization of helical content, active-site disruption and increases in the flexibility of the occluding loop. All these events are closely related to loss of enzyme activity. Heparin binding counteracted these effects by a long-distance, allosteric mechanism which appears to be associated with (i) conformational changes leading to a more compact interface, and (ii) the restriction of catB flexibility, allowing stabilization and rearrangement of interdomain contacts under neutral/alkaline conditions. These new findings may be crucial for the design of new catB inhibitors.
The three-dimensional atomic coordinates of heparin were taken from the Protein Data Bank (entry 1HPN) . From this structure we considered a heparin disaccharide for our studies. For docking and molecular dynamics we used the partial charges and force field parameters for heparin described in ref. . The human cathepsin B crystal structure (PDB entry: 1HUC)  contains two proteins per asymmetric unit; the protein consisting of chains A and C was used as the starting structure for docking and MD simulations. The coordinates of catB complexed with CA-030 (PDB entry: 1CSB)  were also used to increase conformational sampling. We removed the coordinates of the inhibitor and the enzyme was treated as an apo form.
To obtain the pKa values for titratable catB residues we applied the PROPKA web server , using the catB crystal structure as input. This tool uses empirical parameters to estimate pKa values of each titratable residue in its chemical micro-environment, considering the contribution of H-bond formation, charge interactions and desolvation effects. This method was recently compared to other pKa predictive tools obtaining excellent results .
With the results of the pKa analysis, we determine the most probable protonation state for each pH condition and added hydrogen atoms accordingly. Electrostatic surface calculations were performed using APBS software . This analysis combines the solvent accessible surface (SAS) calculation with the values of the electrostatic potential for each atom at the surface. This latter calculation is obtained by the resolution of the linearized Poisson-Boltzmann equation to obtain the surface charge distribution. Images were generated using PyMol .
The heparin disaccharide was docked as a ligand with free rotation for all its torsion angles (sulfate, hydroxyl groups and glycosidic bonds). Docking calculations were performed using the Lamarckian Genetic Algorithm (LGA) implemented in Autodock 3.0 , a program extensively used to predict proteins/polysaccharide binding modes . We applied a two-step protocol. Initially we performed blind docking with a grid of probe-atom interaction energies and the electrostatic potential generated covering the whole protein with a spacing of 0.498 Å. This step was necessary to obtain the putative heparin binding sites on catB structure. Subsequently, we carried out runs taking into consideration only the regions around the most populated low energy clusters obtained from the blind docking. In this step we used a grid spacing of 0.202 Å to evaluate more precisely the atomic interactions in the binding sites. In all, we performed a total of 100 docking runs using a population of 200 individuals. The solutions within 2 Å root mean square deviation of each other were grouped in the same cluster, which were then ranked according to their docking energy values. We selected docked structures of heparin from the lowest energy clusters obtained in our calculations to assemble two distinct catB-heparin complexes, in which heparin was bound to either one or the other of the two catB domains. We identify these complexes according to the domain occupied by the sugar (L or R) and used these structures as starting points for the energy minimization/MD protocols.
Molecular dynamics (MD) simulations, energy minimization and trajectory analyses were carried out with the GROMACS 3.1 package [51, 52] using the GROMOS96 (G53a6) forcefield . Different pH conditions were mimicked by changing the protonation states of six key residues accordingly to pKa calculations and experimental data as described above. Explicit SPC water molecules  were used in all simulations, in which a 14 Å layer of water molecules were added around the solute molecules within a cubic water box, using periodic boundary conditions. Counter ions were inserted for system neutralization. LINCS  and SETTLE  were applied to constrain solute and solvent bonds, respectively. Temperature and pressure were kept at 310 K and 1 atm respectively by the Berendsen approach . Electrostatic interactions were calculated with the PME method , using non-bonded cutoffs of 1.0 nm for Coulomb and 1.2 nm for van der Waals. The MD integration timestep was 2fs.
A three-step energy minimization protocol was used to avoid artifacts in atomic trajectories due to conversion of potential into kinetic energies: firstly, we applied the steepest-descent algorithm: i) 5000 steps with solute heavy atom positions restrained to their initial positions using a harmonic constant of 1 kJ/mol.nm in each Cartesian direction, allowing unrestrained water and hydrogen movement; and ii). 5000 steps with all atoms free to move. Subsequently, the conjugate-gradients algorithm was applied for further energy minimization until an energy gradient fell to 42 KJ/mol.nm. A preliminary MD (1 ns) with heavy atom positions restrained was performed to achieve solvent equilibration and system heating to 310 K. In this step, initial velocities were generated for each simulation. A replica simulation was performed for each system using different initial velocities. Then we performed a 3 ns unrestrained equilibration MD for each system followed by a 40ns production MD.
In the simulations of the heparin-bound complexes, we verified the heparin-catB interaction stability by measuring the time-evolution of the distances between the centers of masses of the heparin and catB positive domains (Additional file 2). This control is important to guarantee that the effects observed result from the maintenance of a stable complex during the entire simulation.
We used the last 20 ns of the catB trajectories to obtain the covariance matrix of C-alpha atomic positions. In this step we applied the g_covar module of GROMACS package. Rotation and translation motions were removed prior to covariance matrix calculation by least-squares superposition. All analyses were performed with the g_anaeig module of GROMACS.
Domain motions were analyzed with the Domain Select module of the DynDom software. This program compares two available structures of the same protein and deduces the rigid-body movement of one domain (the moving domain) relative to the other domain (the fixed domain) in the same way as the DynDom program. We selected the starting structure of each MD and the average structure of the respective simulation to be compared.
We applied the g_cluster module of GROMACS package to calculate the RMS clusters of catB backbone conformations. Our clustering criteria was 1Å RMSD.
We are grateful to MSc Roberta S Faccion for critical reading of the manuscript. We thank FAPERJ, CNPq and CAPES for financial support.
This article has been published as part of BMC Genomics Volume 11 Supplement 5, 2010: Proceedings of the 5th International Conference of the Brazilian Association for Bioinformatics and Computational Biology. The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2164/11?issue=S5.
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.