Structural analysis and insight into novel MMP-13 inhibitors from natural chemiome as disease-modifying osteoarthritis drugs

Purpose: To identify natural chemiome that inhibits matrix-metalloproteinases (MMPs) with a view to discovering novel disease-modifying osteoarthritis drugs (DMOADs). Methods: Computer-aided drug design (CADD) with virtual screening, ADME/Tox, molecular docking, molecular dynamics simulation, and MM-PBSA calculations were used in search of novel natural compounds that inhibit MMPs. Results: From more than fifty thousand compounds, a single lead compound (IBS ID: 77312) was shortlisted using bias based on binding energy and drug-likeness. This lead compound synergistically bound to the S1 domain of MMP-13 protein through five hydrogen bonds. The interactions became stable within 100-nanosecond molecular dynamics simulation run. The in vitro data for the lead compound showed that its minimal non-lethal dose increased collagen content but decreased aggrecan level in chondrocytes. Conclusion: This study has identified a natural lead compound that may pave the way for a novel DMOAD of natural origin against OA.


Osteoarthritis
(OA) is an autoimmune degenerative disorder and the main cause of physical disability and poor quality of life in industrialized countries [1].The disease which affects the joints, is characterized by irreversible damage to articular cartilage, and several studies have so far focused on the molecular biology of the events that lead to cartilage erosion [2].The loss of articular cartilage in OA is triggered by gradual degradation and fibrillation of the cartilage surface.The etiology of articular cartilage degradation lies in the damage to type II collagen [3].Type II collagen, a structural protein which is affected in OA patients, is considered being as one of the main prime targets for new anti-OA drug discovery.There is usually a strict balance between collagen type II production and its degradation by catabolic enzymes.This balance is disrupted due to increased protein degradation [4].Matrix metalloproteinases (MMPs) are catabolic enzymes responsible for collagen type II degradation [5].
Disease-modifying OA drugs (DMOADs) are compounds capable of delaying or halting the progressive destruction of joint tissuein OA.These drugs can be developed through advancements in the understanding of the pathophysiology of OA, and identification of promising therapeutic targets [6].Human-3 collagen (MMP-13) has been widely reported as a target for development of DMOADs [7].In mice deficient in MMP-13, OA cartilage erosion is reduced [8].Human-3 collagen (MMP-13) belongs to the family of zinc-dependent enzymes that are assembled in synovial cells and chondrocytes.It is responsible for decomposing cellular matrix collagen type II protein and maintaining cartilage-specific matrix phenotype [9][10][11].
So far, a few compounds have been reported as DMOADs, but no studies have investigated natural compounds as DMOADs against MMP-13 [12, ].In order to identify DMOADs of natural origin, a natural product database (interBioScreen) that includes about sixty-six thousand compounds was selected.In this database, 60 -65 % of the compounds are of plant origin, 5 -10 % are of microbial origin, 5 % are from marine species, while the rest are from other natural sources.In this study, in silico tools were used to search for natural DMOADs against MMP-13.Virtual Screening (VS) of sixty-six thousand compounds gave < 1 % of the compounds for further evaluation on druglikeness parameters.Eight compounds were subjected to molecular docking simulations against the active site of MMP-13.The top lead compounds were evaluated at atomic level using 100 nanosecond (ns) molecular dynamics simulation with GROMACS atomic simulation tool.The trajectory from the simulations was used for molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA) analysis, for understanding the binding mode of the top lead compound against MMP-13.The top lead compound was exposed to rat chondrocytes under in vitro conditions.

MMP-13
protein in complex with 2napthylsulfonamide hydroxamic acid inhibitor downloaded from RCSB protein databank (PDB ID 3ZXH, 1.30Å x-ray resolution) was used in this study.The structure is a dimer and its monomer from104-272 amino acids was used.The structure was stripped of 2-napthylsulfonamide hydroxamic acid inhibitor and water molecules, and the missing atoms in some of the amino acids in the selected crystallographic structure were added.Energy minimization was done with SPDB viewer.Using OPLS 2001 force field, a final root mean square deviation (RMSD) value of 0.3 Å was obtained.The natural compound database used in this study was downloaded from interBioScreen website (https://www.ibscreen.com/).

Virtual screening for prediction of drug likeliness
Virtual screening was done using ArgusLab suit which limited the number of compounds based on binding energy bias.The S1 binding region of the catalytic domain was selected for screening of natural compounds.For shortlisting, only compounds showing binding energies less than -8 kcal/mol were selected for further studies.Less than 1 % of the total compounds were shortlisted, based on their binding energy (∆G) calculations.The selected compounds were further limited by Lipinki Rule of five (RO5) [13].The RO5 parameters gave eight compounds for molecular docking.

Molecular docking simulations
The shortlisted eight compounds were subjected to atomic scale target-based drug discovery through molecular docking simulation performed using AutoDocksuit (V.4.2.6)[14].AutoDock suit calculates binding free energy for the best protein-compound complex.The suit uses the binding location definition in the form of a grid to create a rectangular network, and the dockings of the compound are usually performed within this rectangular network.AutoDock suit calculates energy values which represent the total energy of the protein, the total energy of the compound and the total torsional free energy.In this study, the results were generated by setting the AutoDock suit to the default parameters, where the maximum ratings were evaluated to illustrate the top ten binding styles per run.Docking results were analyzed using Discovery Studio Visualizer.
The best combination of natural compound with MMP-13 was selected for further assessment using molecular dynamic simulation tool.

Molecular dynamics simulation
The top lead compound for DMOAD against MMP-13 was studied for stability and dynamic interaction using GROMACS atomic simulation platform over 100 nanosecond (ns) time period [15].The protein was subjected to GROMOS 43a1 force field, and force field of the compound was calculated by using online PRODRG web server [16,17].A water bath model with simple point charge (SPC216) was used for water molecules in this simulation.The whole system was energy-minimized by NVT and NPT ensembles, which are governed by Berendsen coupling scheme.For NVT ensemble, a constant temperature of 300 K with 1ns time interval and 0.1 ps coupling constant was used.This step was followed by a NPT ensemble where 1 ns time interval and constant pressure of 1 bar with constant coupling of 5 ps was used.A 100-ns simulation of the complex at 300 K was carried out using SHAKE algorithm and particle mesh Ewald method.In Ewald method, a 14-Å cutoff for van der Waals interactions, and a 12-Å cutoff for Coulomb interaction were updated after every 2 steps [18].Pymol, Ligplus and VMD tools were used for data analysis.The graphs were plotted using Grace program.

MM-PBSA calculations
The MMP-13 and lead complex interaction at residue level were used for MM-PBSA calculations, based on the trajectory from the simulation [19].For MMP-13 and lead complex trajectory, 1000 snapshots were taken, and each snapshot was retrieved after 100ps.The binding mode of MMP-13 with the natural lead compound was studied using grace plotting software.The binding energy (∆G bind ) of MMP-13 protein with the shortlisted natural compound in water was calculated with MM-PBSA using Eq 1.
Where G c is total free energy of MMP-13 and lead complex, while G P and G L are the total free energies of MMP-13 protein and compound in water, respectively.The free energy for MMP-13 and natural compound was obtained as in Eq 2.
where G x represents free energy of MMP-13 or natural compound, E mm is their mean molecular mechanics potential energy in vacuum, TS represents entropic contribution to E mm in vacuum, and G solv is the free energy of solvation.

Isolation and culture of chondrocytes
Female Sprague-Dawley rats weighing about 120 g were purchased from the Laboratory Animal Center of Xi'an Jiaotong University (Xi'an, China).Fetal bovine serum (FBS) and DMEM were purchased from Invitrogen.Hyaluronidase, collagenase, trypsin, and alcian blue 8GX were products of Sigma.Mouse anti-rat collagen monoclonal antibody from Neomarker.Goat antirat aggrecan monoclonal antibody was obtained from Santa Cruz Biotechnology.DAB kit was purchased from Beijing Zhongshan, RNA isolation kit was product of Shanghai Feijie, while first strand CDNA synthetic kit was obtained from MBI.
The isolation and culture of rat costal chondrocytes were performed as described previously [20].The rats were sacrificed by cervical dislocation.Their ribs were removed and the cartilage was collected and cut into 1-mm 3 portions.The tissue samples were digested with hyaluronidase and collagenase, and the isolated cells were cultured at 37°C under a humidified atmosphere containing 5% CO 2 .The cells were sub-cultured at a ratio of 1:1 after treatment with 0.25% trypsin/EDTA.The isolated chondrocytes were identified by immunostaining for collagen and aggrecan.In this study, the cells were used at passage two.

Treatment of chondrocytes with top lead compound
Cells (5 × 10 4 ) were seeded in triplicates in 24well plates and incubated with different concentrations of lead compound (1, 10 and 100 µM) for 48 h.The cells were then collected and tested for expression.

Immunocytochemical analysis
Cells were placed on coverslips and were grown to 80% confluence.They were washed and fixed in 4 % paraformaldehyde for 30 min, and nonspecific binding sites were blocked with normal goat serum.The cells were then incubated at 4°C overnight with anti-aggrecan (1:200 dilution) and anti-collagen (1:200 dilution) antibody, followed by incubation with biotinylated IgG (1:400 dilution).The bound antibody was visualized using 3,3'-diaminobenzidine.Coverslips were mounted and the cells were examined under a microscope.

Drug likeliness
The methods used in this study are shown in Figure 1.The S1 pocket of the MMP-13 protein is the active cleft (site) used for virtual screening of the sixty-six thousand compounds of natural origin.The S1 pocket has the active site cleft which is the continuation of the catalytic domain in β-sheet conformation, forming hydrogen bonds with S1 pocket wall Pro242.The virtual screening grid (Figure 2 b) was created around the S1 pocket and flexible docking was used for shortlisting the compounds.Binding energy score was used to shortlist compounds for further DMOAD development.For this, a binding energy bias was set where only compounds with ∆G less than -9 kcal/mol were shortlisted.The maximum binding energy shown in this screening was -13.765 kcal/mol, and only 582 natural compounds showed binding energies in the range of -9 to -13.765 kcal/mol.The selected compounds were further shortlisted by checking them with parameters set by Lipinski, where four properties i.e. molecular weight (MW), octanol/water partition coefficient (clogP), hydrogen bond donor (HBD) and hydrogen bond (HBA) were evaluated.The method is an established procedure in lead development.Indeed, 90%of the oral drugs that passed phase II clinical trial follow RO5 parameters [17].
The shortlisted compounds were further checked for toxicity, and only compounds which were noncarcinogenic and non-mutagenic were shortlisted.After drug likeliness prediction analysis, eight compounds were subjected to molecular docking simulation analysis (Figure 3).

Molecular docking simulation
Details of the eight compounds selected for atomic scale DMOAD development against MMP-13 protein are given in Table 1.
Each compound was subjected to docking by AutoDock docking tool into the optimized S1 active pocket of MMP-13 protein.Maximum evaluations were performed for each docking and top ten binding poses were generated for each compound.The top generated pose of each compound was studied using Discovery Studio Visualizer and the results are shown in Table 2. Two compounds out of eight formed hydrogen bonds with S1 pocket of MMP-13 protein.

MMP-13 S1 pocket jammers
Two natural compounds (Table 2) were found to jam the S1 pocket of MMP-13 protein with noncovalent hydrogen bond interactions.
The first shortlisted compound was IBS ID: 64787, and its interaction with S1 pocket is displayed in Figure 4.This compound showed ΔG of -8.23 kcal/mol, and it jammed the S1 pocket through formation of two hydrogen bonds with Pro242and Tyr244.The binding pocket for the compound was formed by the amino acids His187, His226, Glu223, Ala186, Leu184, Leu185, His232, His222, Tyr244, Val219, Phe 241, Ile243 and Pro242.The oxygen at first position (O-1) of the compound formed hydrogen bond with Pro242 and the hydrogen bond length was 2.53 A˚.
The oxygen at third position (O-3) of the compound formed one hydrogen bond of 1.92 A˚ length with Tyr244.The mapping of IBS ID: 64787 against the hit (Pro242and Tyr244) is shown in Figure 5a.The second best compound of natural origin was IBS ID:77312 which interacted with the active site of MMP-13 through 5 hydrogen bonds.Its ΔGwas-9.67 kcal/mol, and its binding pocket comprisedHis187, His226, Glu223, Ala186, Leu184, Leu185, His232, His222, Tyr244, Val219, Phe 241, Ile243 and Pro242.The details are shown in Table 3. Figure 5 shows the mapping of the five hits by IBS ID:77312.Based on the molecular docking simulation analysis, IBS ID:77312 was selected as the lead compound for the development of DMOAD.

Molecular dynamics simulation
Molecular dynamics simulation was carried out for the dynamic analysis of the complex to monitorlead-MMP13 complex interaction over time, and 100 ns MDS was done using GROMOS 43a1 force field.The lead-MMP13 complex was stable throughout the run.Root mean square deviation (RMSD) was used to evaluate the MMP-13 protein structure convergence to equilibrium after the interaction of the compound.The backbone atoms of MMP-13 were used to evaluate the RMSD plot (Figure 6a).The complex converged with RMSD range almost above 0.4 nm.The hydrogen bond plot results showed a range of 1 to 7hydrogen bonds formed in 100 ns simulation (Figure 6b).The hydrogen bond formation showed strong binding of the shortlisted compound with MMP-13.

MM/PBSA
The MDS trajectory of the lead compound and MMP-13 protein was used for MM-PBSA to evaluate binding free energy over a period of 100ns.The complex of protein and the compound was taken after every 100-ps of stable intervals from the trajectory and were used to calculate the binding free energy.The MM-PBSA results indicated that the lead compound possessed a negative binding free energy of -244.764kJ/mol.The relative binding free energies of the complex of lead-MMP13 showed strong binding in the dynamic system.The interaction with S1 pocket of MMP13 is shown in Figure 7.

Morphological and immunocytochemical findings
The isolated cells attached to the culture plate within 24 h and showed a typical polygon shape.
After culturing for 3 -5 days, cells grew to confluence and could be sub-cultured.The amount of aggrecan and collagen II (Figure 8) expressions due to increases in the concentration of lead compound showed that the lead compound is a potential candidate for DMOAD development.

CONCLUSION
The shortlisted lead compound can synergistically inhibit the S1 domain of MMP-13 protein.This study has successfully identified a natural compound that may pave the way for a novel DMOAD of natural origin against osteoarthritis.

Figure 1 :
Figure 1: Methodology used to find a novel lead of natural origin against osteoarthritis Figure 2 a shows the pocket selected for virtual screening with a possible DMOAD in the cavity.The virtual screening grid (Figure 2 b) was created around the S1 pocket and flexible docking was used for shortlisting the compounds.Binding energy score was used to shortlist compounds for further DMOAD development.For this, a binding energy bias was set where only compounds with ∆G less than -9 kcal/mol were shortlisted.The maximum binding energy shown in this screening was -13.765 kcal/mol, and only 582 natural compounds showed binding energies in the range of -9 to -13.765 kcal/mol.The selected compounds were further shortlisted by checking them with parameters set by Lipinski, where four properties i.e. molecular weight (MW), octanol/water partition coefficient (clogP), hydrogen bond donor (HBD) and hydrogen bond (HBA) were evaluated.The method is an established procedure in lead development.Indeed, 90%of the oral drugs that passed phase II clinical trial follow RO5 parameters [17].

Figure 2 :
Figure 2: (a) The S1 pocket of MMP-13 protein used to design DMOAD of natural origin; (b) the grid around Si domain used for virtual screening

Figure 4 :
Figure 4: Non-covalent interactions formed by top lead compound IBS ID:77312

Figure 6 :
Figure 6: (a) RMSD plot of the top complex over 100ns simulation, black represents the protein and red represents lead compound IBS ID:77312; (b) The hydrogen bond plot between ligand and protein showing the number varying from one to seven hydrogen bonds over a period of 100ns

Figure 7 :
Figure 7: Binding energy distribution per residue showing the binding of ligand with S1 pocket of MMP-13 protein

Figure 8 :
Figure 8: Effect of different concentrations of lead compound on the expression of aggrecan and collagen II in chondrocytes DISCUSSION OA has unknown pathophysiology affecting joints and leading to physical disability.In-silico techniques have been used to look for drugs for auto immune disorders like rheumatoid arthritis and OA (21,22).Natural compounds have been treated as a goldmine to look for drugs for arthritis {23].OA currently lack treatment options [6] and looking for new strategies and compounds is the need of the hour.Xu et al showed compounds in Curcuma longa, Sophora japonica and Camellia sinensis modulate proinflammatory cytokines TNF-α and IL-1.Khan et al 2017 showed that Wogonin, a compound of natural origin modulates inflammation and shows chondroprotective effects.Chemical database (Maybridge) of 59000 compounds have been used for pharmacophore-based virtual screening to identify novel and potent ADAMTS-4 inhibitors, this study with its limitation to only docking-based studies have reported five possible ADAMTS-4 inhibitors (24).The use of natural compounds with both in-silico and in-vitro validation gives this study an advantage for further biological activity test.

Table 1 :
Drug-likeliness of top Eight lead compounds

Table 3 :
Binding Pocket

Table 2 :
Binding pose analysis of top eight compounds, the ligand binding pocket and the hydrogen bond formation was calculated using Discovery Studio v 16.1.015350