MOLECULAR DYNAMICS STUDIES OF FULL HUMAN MATRIX METALLOPROTEINASE 9 LIGANDED WITH N-HYDROXY-2-[(4-PHENYLPHENYL) SULFONYL-PROPAN-2-YLOXYAMINO] ACETAMIDE

The research presented in this article aimed to provide a full quarternary structure of human matrix metalloproteinase 9 (MMP9) enzyme with a ligand in the catalytic site for structure-based virtual screening. The enzyme plays an important role in wound healing of diabetic foot ulcer. By employing the primary structure of the enzyme obtained from UniProt database (UniProt:P14780), the theoretical structure of full apoenzyme of the human MMP9 (PDB:1LKG), the crystal structures of the catalytic domain (PDB:4H3X) and the hemopexin domain (PDB:1ITV) of the human MMP9, homology modeling studies have been performed. The ligand N-2-(biphenyl-4-yl-sulfonyl)-N-2-(isopropyloxy)-acetohydroxamic acid (CC27) or Nhydroxy-2-[(4-phenylphenyl)sulfonyl-propan-2-yloxyamino]acetamide (IUPAC version) from PDB:4H3X was embedded in the catalytic site of the enzyme. The modeling made use of the modules of homology modeling in YASARA structure. Subsequently, molecular dynamics (MD) simulations in YASARA structure were performed to examine the stability of the enzyme. The homology model was found stable after 5.05 ns and the lowest energy of the model was found at the 6.40 ns of the MD production run. This lowest energy snapshot was then energetically minimized and analyzed for its applicability for virtual screening. This optimized model was then stored in Mendeley Data (DOI: 10.17632/4gsb4p75gz.1).


INTRODUCTION
Enzyme matrix metalloproteinase 9 (MMP9) becomes a molecular target of interest in the discovery of therapeutic agents for a diabetic foot ulcers and cancer (Hariono et al., 2018). The enzyme comprises a catalytic domain and a hemopexin-like domain, which have different roles (Roeb et al., 2002). The catalytic domain degrades the damaged matrix membrane in the wound healing processes (Vandooren et al., 2017). Uncontrolled enzyme activity causes the membrane degradation and formation balance disruption in the wound healing processes (Ayuk et al., 2016). Inhibitors in the catalytic domain are required to control the balance (Jones et al., 2019). On the other hand, although existed in the same MMP9 enzyme, the hemopexin domain has a different role (Dufour et al., 2010). Hemopexin domain acts as a messenger receptor by forming dimers (Ezhilarasan et al., 2009). The interaction between the hemopexin domain with cell surface receptors triggers cell proliferation (Delozier et al., 2011). The growth of cancer cells could, therefore, be inhibited by inhibitors for the hemopexin domain (Alford et al., 2017).
The available 3D structures can assist the visualization of the ligand and protein 2 Roy Gunawan Wicaksono et al.
interactions in structure-based drug design (SBDD) approaches (Wang et al., 2018). The homology comparative modeling method is an approach to predict the 3D structure of protein-based on its amino acid chain. Threedimensional structures in this approach are made using other similar proteins that the 3D structures are already known as the templates (Cavasotto and Phatak, 2009). The processes of building the structure of the 3D protein model include: (1) Other similar proteins with known the 3D structures are used as templates; (2) The amino acid sequences of protein target are aligned to the templates; (3) The 3D structures of the target are constructed based on the 3D structure of the templates; (4) The MMP9 model resulted are then evaluated, validated, and remodelled until appropriate models are obtained (Mart et al., 2000).
There is one publicly available of full MMP9 model (PDB:1LKG). Unfortunately, there is no ligand in the model, which then create difficulties in the binding pocket identification. The aim of the research presented in this article was to provide validated virtual target to discover inhibitors for MMP9. Therefore, homology modeling approaches to construct a 3D structure of full MMP9 with ligands followed by molecular dynamics simulations to examine the stability of the model were performed. The optimized model resulted in this research could be downloaded from Mendeley Data (DOI: 10.17632/4gsb4p75gz.1).

Materials and Instrumentations
The human MMP9 primary sequences (UniProt:P14780) were obtained from The Universal Protein Resource (UniProt; https://www.uniprot.org/). Three-dimensional structures from The Protein Data Bank (PDB; https://www.rcsb.org/) with identity PDB:1LKG, PDB:4H3X, and PDB:1ITV were used as the 3D structure templates of human MMP9. The software mainly used was YASARA Structure version 19.1.27 (license number of 394125786). Computational studies for homology modeling were performed in a workstation with Intel® Pentium® Silver N5000 as the processor, 4 GB random access memory (RAM), and Windows 10 Home 64bit as the operating system. The molecular dynamics simulations and re-docking simulations were performed in a workstation with Intel® CoreTM i5-7500 as the processor, 8 GB RAM, and Windows 10 Professional 64bit as the operating system.

Procedures
The homology model of the human MMP9 was constructed by employing P14780.fasta obtained from UniProt as the source of the primary sequences and PDB:1LKG (sequences 1 to 93 and 445 to 511), PDB:4H3X (sequences 94 to 444), and PDB:1ITV (sequences 512 to 707) from PDB as the templates. The module homology modeling in YASARA Structure was employed to construct the initial full model of the human MMP9. Subsequently the ligands from PDB:4H3X, and PDB:1ITV were incorporated by aligning the chain A of the structures to the initial model and followed by removing the protein parts of the crystal structures. Focused on the side chain residues of the model, the quarternary structure was then subjected to molecular dynamics simulations with fixed atoms in the main chain of residues number 94 to 444 and 512 to 707. All ligands were also fixed during molecular dynamic simulations. The molecular dynamic simulations were performed in YASARA Structure by employing a modified macro from the default macro md_run.mcr from YASARA Structure (Krieger and Vriend, 2009). The modification was done in line 67, from duration='forever' to duration=1000. The molecular dynamics simulations were followed by trajectory analysis using md_analyse.mcr. The final scene file from molecular dynamic simulations resulted from the analysis was unfix and then subjected to energy minimization using energy minimization module. The scene resulted from energy minimization was then saved as mmp9_final.sce. Objects 2 and 3 in the scene were then joined to object 1 by using module Join in YASARA. The object 1 was then 3 stored in YASARA object and pdb formats as mmp9_final.yob and mmp9_final.pdb, respectively. The model of the human MMP9 was then analyzed by employing module check in YASARA structure and AMBER14 as the force field.
The file mmp9_final.pdb was subsequently employed as the starting point for 20 ns MD simulations using YASARA Structure (Krieger and Vriend, 2009) . The default macro md_run.mcr from YASARA Structure was used with modification in line 67, from duration='forever' to duration=20000 and in line 109, from saveinterval=100000 to saveinterval=10000. The MD simulations were followed by trajectory analysis using md_analyse.mcr. The structure of the complex enzyme-inhibitor is considered stable if the deviation of the rootmean-squared-deviation (RMSD) values of at least 5 ns duration of the MD simulations is less than or equal to 1 Å (Liu et al., 2017). The first 5 ns of the MD simulations were considered as the production run. The free energy of binding (G) of all snapshots during the production run was calculated by employing a method published by Prasasty and Istyastono (2019). The snapshot with the lowest G value was then minimized and selected to be analyzed for its applicability for being the target in structure-based virtual screening campaigns by performing 1000 times redock simulations using VINA embedded in YASARA Structure. The docking configuration was embedded to the macro file dock_run_1000.mcr developed in this research. The file could be obtained from Mendeley Data ((DOI: 10.17632/4gsb4p75gz.1). The model was considered as applicable for virtual screening if the RMSD values of at least 95% redocking result of the ligand are less or equal to 2 Å.

RESULTS AND DISCUSSION
The research presented in this article aimed to publicly provide a full quarternary structure of human MMP9 with a relevant ligand in the catalytic site to be used further in drug discovery for diabetic wound healing. Adherence to naming conventions 0.000 > 0 is bad 4 Normality of bond lengths 0.700 < -2 is poor, < -4 is bad 5 Normality of bond angles 0.014 < -2 is poor, < -4 is bad 6 Normality of dihedral angles 0.057 < -2 is poor, < -4 is bad 7 Normality of planar groups 0.573 < -2 is poor, < -4 is bad *) The criteria of "Normality of water locations" was not checked since the model does not contain any water molecule.
4 Figure 1. Visualization of the catalytic site of the homology model (A) and PDB:4H3X (B). The figures were built using YASARA Structure in default mode without shadow. The enzyme, the amino acid residues, the ligands, and zinc atom are depicted as ribbon style, stick style, ball-and-stick style and ball, respectively. Only important residues are shown for the shake of clarity.
The analysis of the human MMP9 homology model showed an acceptable deviation from the average structure in all aspects. The results and the details of the analyzed aspects are presented in Table 1. Based on the results, the model of the tertiary structure of the human MMP9 resulted from the homology modeling is acceptable for further application. However, visual inspection of the interactions between the protein, the ligands, and the cofactors should be performed and compared to those from the reference structures, i.e., PDB:4H3X and PDB:1ITV.
The catalytic site of the model ( Figure  1A) is slightly different from the catalytic site of the crytal structure PDB:4H3X ( Figure 1B). The histidine triad in the model is slightly away from the zinc atom compared to the crystal structure, while the glutamate residue remains similar. Nevertheless, the quantitave approach by aligning the backbone of the model to the backbone of the crystal structure showed the RMSD value of those residues in the model compared to the crystal structure was 0.817 Å. Figure 1 and the alignment show that the model could reproduce the essential interactions for ligands and cofactors with the protein. Unfortunately, the ligands in PDB:1ITV, which was used as one of the templates, are only sulfate ions. Therefore, further investigation in this hemopexin-like domain should be performed.
The reliability of the model to be employed in a further virtual screening campaign will significantly be increased by employing a model with stable protein-ligand interactions (Liu et al., 2017). Therefore, further MD simulations for 20 ns with a snapshot in every 10 ps were performed. Figures 2A and 2B show the RMSD values of the backbone atoms during the simulations and the deviation of these RMSD values in the duration of 5 ns started from the initial point, respectively. As depicted in Figure 2, the MMP9-CC27 complex in this research was considered stable starting from 5.05 ns of the MD simulations. The G values of snapshots in 5 ns after the starting stable point were then calculated to identify the snapshot with the lowest G value. Figure 3 shows that the lowest snapshot was identified at the 6.40 ns of the MD production run.  The MMP9-CC27 complex from the 6.40 ns of the MD production run was then minimized and exported as a pdb file (mmp9.pdb). The applicability checks of this complex for virtual screening by employing 1000 redocking simulations resulted in RMSD values of less than 2.0 Å for all redocking of CC27 poses. The highest RMSD value was 1.274 Å. Thus, the MMP9-CC27 complex (mmp9.pdb) resulted in this research is highly suggested to be used in the further development of structure-based virtual screening protocols to discover drugs targeting human MMP9.