Computational Elucidation of Human β-defensin-2 as a Dual Inhibitor of MMP-9 and PKC-βII for Diabetic Wound Management

1. Introduction

The Diabetic wound (DW) is a late-stage complication of diabetes mellitus (DM) that emanates from interplay amongst inflammation, vasculopathy, neuropathy, and impaired tissue regeneration, all under the background of insulin resistance [1,2]. Although various therapeutic approaches such as pressure off-loading, debridement, and infection control, the management of DW is still in its infancy. The treatment of DW is complex because optimization of glycaemic control has little or no effect; hence, the unique etiopathogenic connotation between DM and DW becomes less relevant at later phases of disease progression [3]. This theoretical problem commands the necessity of developing a novel strategy focused on targeting the ‘specific’ molecular modalities of DW.

The major pathogenic manifestations of impaired healing of DW are an abnormal neutrophil response to the injury and insufficient vascular supply, causing prolonged inflammation and recurrent tissue damage [4]. In general, an extracellular trap (NET) is formed by the active neutrophils with proteins and chromatin granules to protect the wound from infection [5,6,7]. But, in DW, various pathological conditions desensitize neutrophils and favor cellular death by releasing nuclear materials within NETs (NETosis) [8,9].

Gratifyingly, many shreds of evidence exemplified that NETosis and angiogenesis are controlled by the single enzyme protein kinase C-βII (PKC-βII) [10,11,12,13]. Due to excess glycolytic and diacylglycerol production in diabetic conditions, PKC-βII levels are increased and its accumulation results in impaired NETosis and angiogenesis by desensitizing neutrophils and downregulating endothelial nitric oxide synthase, respectively [10,11,12,14,15].

In addition, one more culprit for the prolonged inflammation in DW is a matrix metalloproteinase (MMP)-9, which is highly expressed in inflammatory cells and wound cells [16]. Elevated levels of this enzyme degrade the formed ECM, which delays the DW healing [17,18].

In this scenario, inhibiting PKC-βII and MMP-9 with a single therapeutic agent might be more effective in targeting the multifaceted DW.

Antimicrobial peptides (AMPs) are a part of the innate immune system capable of fighting against infections [19,20]. AMPs have been presented as a potential new antimicrobial class, partly due to their lower susceptibility to bacterial resistance evolution. In addition, they have also been recognized for their potent wound healing activity [21,22]. Human beta defensins (HBDs) are well-categorized AMPs with broad-spectrum activity and are predominantly secreted from leukocytes and epithelial tissues. These are of six types: HBD-1, HBD-2, HBD-3, HBD-4, HBD-5, HBD-6. Among them, HBDs (1-4) are recognized in various parts of the human body, and HBDs (5&6) have been identified recently in human epididymis [23].

HBD-2 is a cationic, antimicrobial peptide discovered in the year 1976, which is a part of innate and adaptive immunity. It comprises 41 amino acids with a gene of 4kb, expressed mainly in the skin, epithelium, and mucous. In general, HBD-2 consists of six cysteine residues at positions 1-5, 2-4, and 3-6, united by three disulphide bridges [24,25,26]. The secondary structure of HBD-2 reveals the presence of an α-helix, three β-sheets organized in an antiparallel mode in N-terminal region and C-terminal region (Figure 1).

The basic folding of the peptide is possible by the disulphide bonds between the α-helix and first β-sheet, which possess a domain of Gly-X-CysIV. Though, the active site present in the N-terminal region, antimicrobial activity of the peptide is mainly because of the presence of cationic residues present in the C-terminal region of the peptide. Remaining biological properties of the peptide is dependent on specific conformation of N-terminal region of the HBD-2. Further, its amphipathic nature and the presence of disulphide bridges is responsible for resistance to degradation by proteases [27,28]. Hence, HBD-2 a biological macromolecule is selected as a ligand for the preliminary computational study with an anticipation to accelerate DW healing by inhibiting PKC-βII and MMP-9 enzymes.

2. Materials and Methods

2.1. Protein-Protein Docking Studies

The The protein 3D X-ray crystal structures of MMP-9 and PKC were retrieved from RCSB-PDB database with PDB IDs of 4XCT and 2I0E. Each receptor structure was validated using PROCHECK programs and ProSA-web [29,30]. Energy refinement was achieved by adding hydrogen atoms and assigning partial charges.

The 3D X-ray crystal structure of HBD-2 was retrieved from the database and considered as ligand in the docking studies. Protein-protein docking was performed by discovery studio using default parameters (2000 protein poses and 15 angular step size) Z-DOCK protocol (http://zdock.umassmed.edu/). Z-DOCK is usually a grid-based algorithm using novel shape-complementarity scoring function. To start the protocol, parallel processing was kept at “false” and Z-Rank at “True”.

2.2. Molecular Dynamics Simulation (MDS) Studies

The complex structures of MMP-9 and PKC-II with HBD-2 as a ligand protein identified from Z-Dock pose were subjected to MDS studies. When compared with docking studies, MDS studies provide better insights of ligand binding [31]. System preparation of both the complexes was carried out using the procedure reported by Gangadharappa et al., 2020 and Prasanth et al., 2021 [32,33]. GROMACS simulation package was utilized to analyze the simulated data of two complexes. Root mean square deviation (RMSD) and root mean square fluctuations (RMSF) were employed to analyze the backbone atoms and C-α fluctuations, respectively. Further, radius of gyration was computed to determine the “extendedness” of the protein residues. All these calculations were performed using various utilities like “gmxrms”, “gmxrmsf”, “gmx select”, and “gmx gyrate.”

2.3. Free Energy Calculation (MM-PBSA)

It is the method employed to obtain insights into the protein-protein interactions by calculating interaction free energies. The free energy was computed using the g_mmpbsa, utility of the gromacs [34]. To validate the MD result, the last 500-1000 ns of the MD trajectories were chosen to obtain free energy.

2.4. Principal Component (PC) and Free Energy Landscape (FEL) Analysis

PC was usually carried out to generate mass-weighted covariance matrix (CM) of the protein atom displacement, indicates the collective and dominant modes of the receptor throughout the simulation studies of the MD trajectory [35,36]. CM is diagonalized to extract a set of eigenvalues and eigenvectors using g_covar (GROMACS utility), which indicates the concerted molecule motion. In addition, eigenvectors were analyzed and plotted using g_anaeig tool, available in GROMACS [37,38]. Further, the results obtained from PC analysis were subjected to FEL analysis. The g_covar and g_anaeig utilities were employed to identify the changes in the motion patterns of both the complexes [39,40].

3. Results and Discussion

3.1. Protein-Protein Docking Studies

Protein-protein docking was performed for HBD-2 (PDB ID: IFD3) against two receptors MMP-9 (PDB ID: 4XCT) and PKC-βII (PDB ID: 2I0E) to recognize the potential ligand for accelerated healing of DW. Although so many conformations were generated (Table 1), the top scored conformation was given preference and used for the further analysis. PDBsum and Pymol visualization tools were used analyse protein-protein interactions. Depending on the binding interactions, two conformations were selected and has exhibited a significant Z-score -118.026 (4XCT) and -140.319 (2I0E).

3.2. Molecular Interaction Studies of Top-Ranked Complex HBD-2/4XCT.

The protein-protein interactions were tabulated in Table 2. The amino acid residues involved in the protein-protein interactions were obtained using the PDB sum of the 3D structure of the complex HBD-2/4XCT. From the analysis, it was clear that only 15 residues of HBD-2 (chain B) and 19 residues of 4XCT (chain A) (Figure 2a). It was observed that the complex exhibits a total of 6 hydrogen-bonding and 237 non-bonded interactions (Figure 2b). Amongst six hydrogen bonds, only two were found to be prominent as depicted in Figure 2c. The remaining hydrogen bonds were not visible as the distance is greater than 2.5 Å. The >C=O and -NH group of Tyr248 exhibited two hydrogen bonding interactions with -NH2 (>C=O⋯HN(H), 1.7 Å) and >C=O (-NH⋯O=C<, 2.2 Å) group of Lys279, respectively. Moreover, the complex was stabilized by the occurrence of non-bonded contacts between the residues.

3.3. Molecular Interaction Studies of Top-Ranked Complex HBD-2/2I0E.

The protein-protein interactions were depicted in Table 3. The protein-protein interacting residues were analysed using the PDBsum of the 3D structure of the top-ranked complex HBD-2/2I0E. During analysis, it was observed a total of 18 residues of HBD-2 (chain B) and 22 residues of 2I0E (chain A) were involved in the interactions (Figure 3a). Further the complex exhibited 2 hydrogen-bonding, 4 salt-bridge and 154 non-bonded contact interactions within the catalytic pocket of 2I0E (Figure 3b). The two hydrogen bonds were found to be prominent and was depicted in Figure 3c. The >C=O of Lys468 (chain A) showed a hydrogen bond with the side chain -NH2 group of Gln695 (chain B) (>C=O⋯HN(H), 1.5 Å). Another hydrogen bond was observed between the >C=O of Gly543 (Chain A) and -NH group of Lys679 (chain B) (-NH⋯O=C<, 2.1 Å). In addition, the complex was found to stabilized due the formation of non-bonded contacts and salt-bridge interactions between the residues.

3.4. Molecular Dynamics (MD) Simulation Analysis of Two Complexes HBD-2/4XCT and HBD-2/2I0E.

A 1000 ns MD simulation studies were performed for each of the docked complexes to explore the protein-protein complex stability and dynamic properties. Average values of every parameter of the two complexes HBD-2/4XCT and HBD-2/2I0E was calculated from their MD trajectories has been tabulated in Table 4. GROMACS simulation package was used to analyse the simulation data. Various gromacs utilities such as “gmxrms”, “gmx select”, “gmrmsf”, and gmx gyrate were used to compute RMSD (backbone atoms), RMSF (Cα) and radius of gyration, respectively.

RMSD images depicts the moving parts and spatial displacement rate of the protein model during MD simulation and also serves as a measure of receptor stability.

The lowers values of RMSD corresponds to protein stability during the simulation and the acceptable average value for a protein is ≤ 0.4 nm. It was observed that average RMSD values for both the complexes were less than 0.4 nm (Table 4). RMSD diagrams were plotted for 1000 ns of MD study with both the complexes HBD-2/4XCT (Figure 4a) and HBD-2/2I0E (Figure 4b). In both the complexes, the systems attained equilibrium after 200 ns. The average RMSD value of complex HBD-2/4XCT was found to 0.303±0.01 nm whereas for the complex HBD-2/2I0E, the value was 0.399±0.06 nm.

Further, RMSF was calculated to determine the dynamic behaviour of the Cα atoms and backbone flexibility of the protein, which in turn useful to investigate each motion that corresponds to the general motion of the protein. The higher values of RMSF correspond to the greater flexibility of the backbone residues during the simulation. RMSFs were plotted for both the complexes HBD-2/4XCT (Figure 5a) and HBD-2/2I0E (Figure 5b) for the 1000 ns MD study. The average RMSF values were found to be 0.338±0.098 nm and 0.451±0.330 nm for HBD-2/4XCT and HBD-2/2I0E, respectively (Table 4).

In case of HBD-2/4XCT, Asp185 to Pro310 lies in the catalytic pocket or protein-protein interface. A total of seven hydrogen bonding interactions were observed between the residues of HBD-2 and 4XCT (Figure 6a). Whereas, in the HBD-2/2I0E, protein-protein interface was occupied with residues Ile671 to Pro710. Approximately eight hydrogen bonding interactions were observed at the interface of MD trajectory pose of the complex HBD-2/2I0E (Figure 6b).

Radius of gyration (ROG) parameter depicts the protein structure compression or folding behaviour during MD simulation. The lower ROG values 1.633±0.034 nm and 2.063±0.024 nm was obtained for the HBD-2/4XCT and HBD-2/2I0E, respectively (Table 4) indicating no conformation of distortion and increased compactness during the MD study (Figure 7a,b).

Further, protein-protein solvent-accessible surface-area (SASA) was computed for both the complexes and represented in Figure 8a (HBD-2/4XCT) and 8b (HBD-2/2I0E). As SASA indicates the solvent replacement phenomenon, lower values indicate that the catalytic pocket is less solvent-exposed and the ligand retains within the catalytic site throughout the study. The average SASA values of 101.278±0.640 nm2 and 183.409±21.349 nm2 was obtained for the HBD-2/4XCT and HBD-2/2I0E, respectively (Table 4), indicating that the HBD-2 has been deeply resided within the binding pockets of 4XCT and 2I0E.

3.5. Binding Free Energy Calculations-MMPBSA Approach (∆Gbind)

The ∆Gbind was computed for both the complexes for each 100 ns of the last 500 ns MD trajectories. The average values of ∆Gbind, Van der Waals (Vdws), electrostatic (E), polar solvation (PS) and SASA energies were computed and tabulated in Table 5. The average ∆Gbind was found to be -2318.73±35.46 kJ/mol for HBD-2/4XCT and -2329.62±61.24 kJ.mol for HBD-2/2I0E complexes, indicating the greater stability of HBD-2 within the catalytic pocket of the both the residues. In addition, the complexes possess significant VdWs, E, SASA and PS, indicating complete burial of HBD-2 and affinity within the catalytic pockets of 4XCT and 2I0E.

3.6. Principal Component (PC) and Free Energy Landscape Analysis

PC analysis is a powerful approach to understanding the motion changes of complexes during MD simulation studies. Diagonalization of the covariance matrix of the Cα atom fluctuations can be analyzed through eigenvalues and was represented in Supp-Figure S1a and S1b, with eigenvector indices. In the case of both the complexes HBD-2/4XCT and HBD-2/2I0E first few eigenvalues relative to the concerted motions, define essential subspace, where the majority of the motions take place. In addition, the eigenvalues were observed to decrease in amplitude in order to reach the more localized fluctuation and constrained numbers (Figure 9a,b).

Further, the deviations found during PC analysis were assessed by computing free energy landscape (FEL) studies. FEL is mainly to evaluate the conformational changes in the biomolecule corresponding with the free energy obtained during the protein-protein interaction and complex formation. The blue regions (higher ∆Gbind) favor an unfolding state of protein whereas the red spots indicate the folding state of the protein. From the FEL plots, it is evident that the two complexes exhibited higher ∆Gbind values or blue spots (Figure 10a,b).

Furthermore, these complexes equilibrated with high coverage of eigenvectors. Almost >95% of total fluctuations were because of the residues in the protein-protein interface. These results can be correlated with the physiological condition as both complexes might be formed more frequently and have attained greater degrees of stability. To visualize a simple, conformational landscape, different frames with a similar conformation were analyzed. The RMSD to native bound structure is also calculated from the 1000 ns trajectory (Supp-Figure 2a and 2b).

4. Conclusions

In summary, the structural properties of interactions, stability and dynamic features of HBD-2 with MMP-9 and PKC β-II were characterized using computational strategies. The results suggest that the amino acid residues present in protein-protein interface which participated in the interactions are vital for the functional activity of the receptors. The docking results were astonishing as the both the top-ranked poses significant Z score. The binding affinity was supported by MD simulation studies. The RMSD, RMSF, SASA, ∆Gbind, PC analysis and FEL analysis clearly demonstrated the structural and behavioral changes of the residues during MD study, which can directly correlate to physiological behavior of proteins. Henceforth, HBD-2 can be directly used as dual-inhibitor of MMP-9 and PKC β-II and developed as a therapeutic for accelerating the DW healing.