Quantum Mechanics/Molecular Mechanics Study on Caspase-2 Recognition by Peptide Inhibitors

For a variety of biological and medical reasons, the ongoing development of humane caspase-2 inhibitors is of vital importance. Herein, a hybrid (Quantum Mechanics/Molecular Mechanics – QM/MM), two-layered molecular model is derived in order to understand better the affinity and specificity of peptide inhibitor interaction with caspase-2. By taking care of both the unique structural features and the catalytic activity of human caspase-2, the critical enzyme residues (E217, R378, N379, T380, and Y420) with the peptide inhibitor are treated at QM level (the Self-Consistent-Charge Density-Functional Tight-Binding method with the Dispersion correction (SCC-DFTB-D)), while the remaining part of the complex is treated at MM level (AMBER force field). The QM/MM binding free energies (BFEs) are well-correlated with the experimental observations and indicate that caspase-2 uniquely prefers a penta-peptide such as VDVAD. The sequence of VDVAD is varied in a systematic fashion by considering the physicochemical properties of every constitutive amino acid and its substituent, and the corresponding BFE with the inhibition constant (Ki) is evaluated. The values of Ki for several caspase-2:peptide complexes are found to be within the experimental range (between 0.01 nM and 1 μM). The affinity order is: VELAD (Ki = 0.081 nM) > VDVAD (Ki = 0.23 nM) > VEIAD (Ki = 0.61 nM) > VEVAD (Ki = 3.7 nM) > VDIAD (Ki = 4.5 nM) etc. An approximate condition needed to be satisfied by the kinetic parameters (the Michaelis constant – KM and the specificity constant – kcat/KM) for competitive inhibition is reported. The estimated values of kcat/KM, being within the experimentally established range (between 10–4 and 10–1 μM–1 s–1), indicate that VELAD and VDVAD are most specific to caspase-2. These two particular peptides are nearly 1.5, 3 and 4 times more specific to the receptor than VEIAD, VEVAD and VDIAD respectively. Additional kinetic threshold, aimed to discriminate tightly bound inhibitors, is given.

The first identified mammalian member is caspase-2 and its physiological role is not quite clear. Caspase-2, one of the most evolutionarily conserved caspases, is inclined to behave as either executioner or initiator. In terms of substrate specificity, caspase-2 is similar to caspase-3 and -7 (executioner caspases). However, the long N-terminal caspase recruitment domain (CARD) of caspase-2 indicates its potential role as an initiator caspase. 3 While the function of caspase-2 in the embryonic development of mice is questionable, 4 its important role in stress-induced cell death pathways and tumor suppression is more certain. 5 The potential roles of caspase-2 in mediating nonapoptotic pathways (cell-cycle regulation and DNA repair) have been reported in terms of whether caspase-2 is mandatory for apoptosis under specific circumstances, 6,7 or whether it primarily functions in cell-cycle regulation. 5 An elevated expression level of caspase-2 has been observed in the brain of patients with some neurodegenerative disorders. 8 In addition, the critical role of caspase-2 in mediating nonalcoholic steatohepatitis (NASH) pathogenesis, a chronic and aggressive liver condition not only in mice but probably in humans, has been highlighted. 9 Whereas many questions on caspase-2 physiology remain enigmatic, one of the key aspects for developing caspase-2 specific probes is related to the way in which caspase-2 gets activated.
Peptide bonds are hydrolyzed using caspases (endoproteases) in a reaction that depends on catalytic cysteine residues in the caspase active site and occurs only after certain aspartic acid residues in the substrate. Besides resulting in substrate inactivation, caspase-mediated processing may generate active signaling molecules that participate in apoptosis and inflammation. Caspase activities are strictly regulated by protein-protein interactions and by proteolysis. 1 The crystal structures of caspase-2 in complex with several peptide inhibitors and comparison of the apo (substrate-free) and inhibited caspase-2 structures have revealed a recognition via several discrete catalytic steps: (i) activation of caspase-2 by breaking a nonconserved salt bridge between Glu217 (caspase-2 is the only human caspase with glutamate at position 217) and the invariant Arg378, (ii) formation of a catalytically competent conformation upon binding to a single substrate, and (iii) formation of the enzyme-substrate complex after having both active sites occupied by the substrate. 10 Caspase-2 has been suggested to uniquely prefer a penta-peptide rather than a tetra-peptide, as required for efficient cleavage by other caspases. 10 To gain more complete insights into the caspase-2/peptide recognition and further facilitate the design of caspase-2 inhibitors, a hybrid QM/MM approach is employed in this work.

Methods
To obtain the initial atomic coordinates of the apo and inhibited caspase-2 structures, the experimental structures were retrieved from the Research Collaboratory for Structural Bioinformatics (RCSB) Protein Data Bank (PDB): 3R7S.PDB (apo caspase-2), 3R6G.PDB (caspase-2/ VDVAD), 3R5J.PDB (caspase-2/ADVAD), 3R7B.PDB (caspase-2/DVAD), and 3R6L.PDB (T380A/VDVAD). 10 The sequence of the penta-peptide inhibitor was varied using single point mutations generated by applying the Mutagenesis engine of PyMol-v0.99 to the experimental structure 3R6G.PDB in a backbone-dependent fashion. 11 Before running QM/MM calculations, the systems were prepared using the Amber 11 suite of programs. 12,13 The solute was prepared using the Amber11 utility program tLeap in association with the ff99sb force field. 14 Every inhibitor was initially prepared by parameteryzing its atom types, charges, and connectivity in order to be treated as part of the solute. The molecular geometry was optimized by Gaussian 98 at the MP2/6-31G* level of theory. 15 The molecular electrostatic potential was calculated by Gaussian 98 at the HF/6-31G* level of theory, 15 while the atomic charges were derived by means of the RESP fitting technique, 16 which is part of AmberTools 1.5. 12,13 Remaining parameters were assigned from the General Amber Force Field (GAFF), 17 being entirely compatible with the ff99sb macromolecular force field. 14 Every solute was solvated using a 10 Å (1 Å = 10 -10 m) pad of TIP3P water molecules (≈ 11500) and the counter ions Na + were added to neutralize each system. To remove clashes and bad contacts, two-stage geometric minimization was performed using the Sander module of Amber11. 12,13 At the outset, the positions of the solute atoms were kept fixed, while the positions of the water atoms were minimized by gradually reducing an initial harmonic restraint of 2 kcal mol -1 Å -2 on all non-hydrogen non-water atoms via 5000 combined steepest descent (2500 steps) and conjugate gradient (2500 steps) minimization steps. Afterwards, the entire system was minimized without restrains by means of 10000 combined steepest descent (5000 steps) and conjugate gradient (5000 steps) minimization steps.
A two-layered hybrid approach was employed to assess the binding affinities within the caspase-2:peptide complexes. The outer layer of the complex (Figure 1) was kept at the low level of theory (MM) with an Amber force field. The central layer of the complex (bold sticks; Figure  1) was treated at the high level of theory (QM), employing SCC-DFTB-D, the Self-Consistent-Charge Density-Functional Tight-Binding method 18,19 with Dispersion energy, 20 as implemented in Amber11. 21,22 The inclusion of the empirical correction for dispersion energy into SCC-DFTB provided a balanced and reliable description of the interactions inside the systems. 23 Pure Density Functional Theory (DFT) methods are known for their modest computational costs, but they are not able to adequately describe dispersive forces, especially within unconventional systems, 24,25 as many density functionals are empirical. 20 DFT was extended to include dispersion correction (DFT-D), 26,27 and as such DFT-D became suitable for performing energy minimizations and vibrational analyses of extended molecular complexes containing hundreds of atoms. By being a few orders of magnitude faster than DFT-D, 23,28-33 SCC-DFTB-D was suggested to be applicable to both quantumchemical simulations and calculations pertaining to a large number of extended molecular complexes. 23 The total interaction energies were defined as: (1) where ΔE model denote the energies of the model system defined at high and low level of theory and ΔE real denotes the whole (real) system. Therefore, the equivalent binding free energies of the complex systems were determined as: The thermodynamic quantities (enthalpies, entropies, and different entropic contributions) were obtained from frequency calculations done by the Nmode module of Amber11. 12,13 The different entropic contributions (translation, rotation, and vibration) for caspase-2:peptide complexes were calculated as:

Results and Discussion
Kinetic measurements of competitive inhibition associated with the initial experimental structures 10 are summarized in Table 1. The specificity constant (k cat /K M ) identified the penta-peptide VDVAD as a preferred inhibitor, while two residues, Thr380 and Tyr420, were identified as critical for recognizing a residue at the P5 position -the first position at the left end in the peptide sequence (Figures 2b & 2c). The salt bridge between Glu217 and Arg378, which is present in the apo caspase-2 (3.37 Å, Figure 2a), is broken in the caspase-2:VDVAD complex (8.05 Å, Figure 2b), because Thr380 and Tyr420 in P5 rec-ognition move 2.1 and 3.6 Å, respectively (Figure 2d). Furthermore, the specificity constant revealed that mutation of Thr380 to Ala reduces the catalytic efficiency of caspase-2 by about 40 fold (Table 1), as Thr380Ala ( Figure  2c) causes the loss of the hydrogen bond between Thr380 and the P5 side chain (3.51 Å, Figure 2b) due to a 2.3 Å movement in the main chain in residue 380. Structurally speaking in a similar manner, mutation of Tyr420 to Ala reduces the catalytic efficiency of caspase-2 by about 4 fold (Table 1), as Tyr420Ala causes a 0.5 Å movement of the side chain of residue 420 and the loss of the hydrophobic interaction between Tyr420 and the P5 side chain. 10 Table 1. Kinetic data 10 for experimental caspase-2:peptide inhibitor complexes: K M -Michaelis constant, k cat -catalytic constant, k cat / K M -specificity constant, and IC 50 -inhibitory concentration  (d) It is unreliable to measure the catalytic efficiency values for the slowest (k cat /K M < 10 -4 μM -1 s -1 ) and fastest (k cat /K M > 10 -1 μM -1 s -1 ) reactions. 36 To perform physically realistic QM/MM calculations, the first important aspect is how to define a QM region, or what caspase-2 residues need to be included in the QM region. There are no good universal rules here. Binding site residues of caspase-2 that are involved in noncovalent interactions with a peptide inhibitor are: Arg219, His277, Gly278, Gln318, Cys320, Ala376, Arg378, Asn379, Thr380, Trp385, Arg417, Glu418, and Tyr420. 10 Caspase-2 is the only human caspase with glutamate at position 217 forming a salt bridge with Arg378 in the apo caspase-2 ( Figure 2a). The inhibition of caspase-2 was related to breaking the Glu217-Arg378 salt bridge, while residues Thr380 and Tyr420 were pointed out as the key elements for recognizing a preferred penta-peptide along a catalytic pathway (Figure 2b). 10 An intention to define a QM region to mimic the active site has to take into account all these experimental and structural arguments. Knowing that inclusion of a different number of caspase-2 residues in the QM region is associated with different thermodynamic   properties such as the binding free energies means that an appropriate QM region is supposed to generate results in agreement with experimental data. Even though one might want to have as large a QM region as possible, having more than 80-100 atoms in a QM region lead to simulations that are computationally very expensive. 12,13 To reconcile all the structural, functional, and computational standpoints as much as possible, the present choice of including Glu217, Arg378, Asn379, Thr380, Tyr420, and the peptide inhibitor in the QM region ( Figure 1) is carefully made in order to generate the inhibition constants that are within an experimental range -from 1 µM (1 µM = 10 -6 M) to 0.01 nM (1 nM = 10 -9 M) for inhibited caspase-2 structures. 34 QM/MM binding free energies for the experimental structures are given in Table 2. Figure 3 shows quite a satisfactory linear correlation between the calculated ΔG bind ( Table 2) and the experimental inhibitory concentration IC 50 (Table 1): ΔG bind = 0.013 IC 50 -12.85, R = 0.97. The most negative BFE for the caspase-2:VDVAD complex (-13.22 kcal mol -1 ) signifies that VDVAD is a favorable inhibitor. The enthalpy contribution (ΔH) for the complexes ranges from -31.85 to -24.28 kcal mol −1 , indicating that the noncovalent complexation process is exothermic. In case of the entropy contribution (TΔS), the less negative entropy change is, the more reduced degrees of freedom of an inhibitor in the protein active pocket are. The least negative entropy (-18.63 kcal mol −1 ) is associated with caspase-2:VDVAD, of which vibrational contribution (4.88 kcal mol −1 ) makes a most conspicuous difference with respect to the other complexes ( Table 2). The increased and thermodynamically favorable vibrational entropy change upon binding of VDVAD to caspase-2 is the signature of preferred noncovalent complexation.
To search for more effective penta-peptides, the sequence of VDVwAD is systematically varied by means of single point mutations of its constitutive residues. To make such a procedure consistent, each amino acid is mutated to its counterpart observed from a physicochemical standpoint. Val (V), an aliphatic and hydrophobic amino acid, is mutated to either Ile (I) or Leu (L). Asp (D), a polar and   negatively charged amino acid, is mutated to Glu (E). Ala (A), a tiny and hydrophobic amino acid, is mutated to Gly (G). The estimated BFE and K i for each generated complex are reported in Table S1 (Supplementary Material). On the basis of the relation ΔG bind = RT ln(K i ) (R -the gas constant = 1.9872 kcal K −1 mol −1 ), T -the absolute temperature = 300 K), the experimental range of K i in between 1 µM and 0.01 nM corresponds to the BFE (ΔG bind ) range in between -8.23 and -15.09 kcal mol -1 for inhibited caspase-2 structures. Numerical inspection of the data in Table S1 identified the complexes that have the BFEs inside of this specific range (Table 3). Thus, as far as affinity issue for the receptor is concerned, the order of preferred inhibitors is: VELAD (K i = 0.081 nM) > VDVAD (K i = 0.23 nM) > VEIAD (K i = 0.61 nM) > VEVAD (K i = 3.7 nM) > VDIAD (K i = 4.5 nM) etc. In order to evaluate the specificity constant for the complexes (Table 3), the correlation of the Michaelis constant with the specificity constant for the experimental caspase-2:peptide structures (Table 1) is observed. Even though two linear correlations are established (Figure 4), the first one (K M = -0.58 k cat /K M + 157.51, R = 0.74; Figure  4, top) is slightly more suitable because it reproduces the experimental value of the specificity constant for the caspase-2:VDVAD complex more accurately than the second one (Figure 4, bottom). Due to the negative slope (-0.58) of the linear regression line, K M < 157.51 µM rep-resents an approximate condition for the physically meaningful estimate of k cat /K M .
To evaluate the functional efficiency of the complexes ( Table 3) in terms of the specificity constant (k cat /K M ), a competitive inhibition mechanism is considered ( Figure 5). For such a reaction, the inhibition constant is defined as: where IC 50 is the inhibitory concentration, K M is the Michaelis constant, and [S] is the substrate concentration. 35 Solving Eq. 4 for the Michaelis constant gives: The kinetic data are analyzed as follows. IC 50 is evaluated using its linear correlation with ΔG bind (Figure 3). K M is evaluated using Eq. 5 with [S] ≈ 2.7 mM (1 mM = 10 -3 M) -a typical experimental value. 10 Of these complexes (Table 3), those having K M < 157.51 µM are selected (Table 4) and may be considered as competitively inhibited structures. The comparison of the values of K M (  (Table 4) is made by way of the linear correlation K M = -0.58 k cat /K M + 157.51 (Figure 4, top). The specificity constants (k cat /K M ) are within the experimental range (between 10 -4 and 10 -1 μM -1 s -1 ), 36 indicating that VELAD and VDVAD are most specific to caspase-2. These two particular peptides are approximately 1.5, 3 and 4 times more specific to the enzyme than VEIAD, VEVAD and VDIAD respectively.
In case of simple enzyme reaction with one substrate, if k cat << k -1 , K M is conceivable as the dissociation constant that quantifies the strength of the ES complex formation. If the K M value gets smaller then the ES complex gets stronger (or more stable). In other words, the more pronounced enzyme affinity to the substrate is lined up with the more specific inhibition of the enzyme. The values of ΔG bind (Table 3), K M (Table 4) and k cat /K M (Table 4) conform to the trend.  For tightly bound inhibitors, the inhibition constant is: (6) where [E] is the enzyme concentration. 35 For tightly bound inhibitors, the equation 6 takes into account the larger amounts of inhibitor bound species, thus making that the Michaelis-Menten assumption of the total enzyme concentration being equal does not hold. 35,37 The Michaelis constant for tightly bound inhibitors is: The condition K M < 157.51 µM for tightly bound peptides means: Solving Eq. 9 for IC 50 gives: that VELAD, VDVAD, VEVAD and VDIAD may be considered as tightly bound inhibitors.
In most experimental investigations of enzyme kinetics, the total concentration of substrate is in excess of the enzyme concentration, thus making the free and total substrate concentrations essentially equal. 35 For a set of inhibitor candidates, comparison of K M or IC 50 values is only assumed to be valid when these values are evaluated under identical experimental conditions. 35,38 Only few data exist on the catalytic efficiencies of caspase substrates, so that a more complete understanding of their true substrate preferences is impossible. 36 The present results are imagined to facilitate additional experiments, which are needed to understand better the kinetics of caspase-2/peptide recognition for further research or therapeutic product development. The fact that caspase inhibition-based drug has not been approved on the market so far means that the development of therapeutic approaches that specifically target caspases is a substantial challenge of particular biological and clinical interest. 39

Conclusions
QM/MM model derived and exploited in this work has been shown to correlate with the existing experimental observations to an appreciable extent, indicating that caspase-2 uniquely prefers a penta-peptide such as VDVAD.
This approach has enabled the extensive and systematic investigations of some of the important aspects both of the thermodynamics and of the kinetics of caspase-2 recognition by a large number of penta-peptides. The sequence of VDVAD has been consistently varied and the corresponding binding free energies with the inhibition constants have been evaluated. The values of the inhibition constants, being within the experimental range for several caspase-2:peptide complexes, have indicated that the affinity order is: VELAD > VDVAD > VEIAD > VEVAD > VDIAD etc. The specificity constants for competitive inhibition have been estimated to fit the experimentally predicted range, thereby suggesting that VELAD and VDVAD are most specific to caspase-2; also both are about 1.5, 3 and 4 times more specific to the receptor than VEIAD, VEVAD and VDIAD respectively. An approximate kinetic   Table 3, those having K M < 157.51 µM are reported here. K M < 157.51 µM is an approximate condition for the physically meaningful estimate of k cat /K M (Figure 4, top).
threshold, supposed to discriminate tightly bound peptide inhibitors, has been reported. This study has demonstrated that a well-calibrated computational work may yield information inaccessible by other methods or suggest new experimental procedures.