A molecular model of human Lysyl Oxidase (LOX) with optimal copper orientation in the catalytic cavity for induced fit docking studies with potential modulators

Lysyl oxidase (LOX) is a copper dependent amine oxidase which catalyses the cross linking of collagen and elastin towards the maturation of extracellular matrix. The expression and activity of LOX is known to vary under pathological conditions such as tumorigenesis, hyperhomocysteinemia, copper deficiency diseases, pseudoexfoliation syndrome and proliferative diabetic retinopathy. Despite the implication of LOX in many diseases, there is inadequate information about its structure. Therefore, we describe a molecular model of Human Lysyl Oxidase (LOX) with optimal copper orientation in the catalytic cavity for induced fit docking studies with potential modulators. The predicted model was found to be highly plausible as per the stereochemistry checks. Further, Molecular Dynamics (MD) studies also inferred the stability of the predicted structure. We performed Induced Fit Docking (IFD) of LOX modulators to the predicted structure and also validated the molecular interactions in implicit solvent model by calculating Molecular Mechanics Generalized Born Surface Area (MMGBSA). The IFD results strongly reveal that aspartic acid residues in the catalytic cavity as the key players in establishing interactions with small molecules. The insights from this study will aid in better exploration of the structure-function relationship of LOX.

LOX is also involved in a spectrum of biological functions which include developmental regulation, tumour suppression, cell migration, adhesion, apoptosis, and cellular senescence [3]. LOX expression is known to be regulated by hypoxia inducible factor-1, transforming growth factor β, tumour necrosis factor α, platelet derived growth factor and fibroblast growth factor [4]. LOX levels are reported to be increased in hepatic fibrosis, Wilson's disease, liver granuloma, lung fibrosis, cardiovascular disease, metastatic breast cancer, plaque deposits in Alzheimer's disease, keloids and scleroderma [1]. Its expression is decreased in Type IX Ehlers Danlos syndrome and Menkes disease [5]. It is well established that LOX is involved in many of the normal biological functions and also in diverse pathophysiological conditions.
The total sequence length of human LOX (Uniport Accession ID: P28300) is 417 residues comprising of three region; A) signal peptide , B) propeptide (22 -168) and C) mature LOX (169 -417). Enzymatically active mature LOX is 249 amino acid residues in length, devoid of signal peptide and propeptide. Copper serves as a prosthetic group in the catalytic domain of LOX. Here, copper exhibits a unique binding pattern when compared with other amine oxidases. In the case of mature LOX, the copper ion is found to be harboured in the talon shaped loop, whereas in other amine oxidases it spans inbetween two beta sheets. The talon loop structure in LOX is comprised of four histidine residues, of which three contribute for coordinated covalent bond interactions with the copper ion The normal and pathological functions of LOX and its isoforms has been the subject of many research over the past two decades. However, there is a significant dearth of information reported on the 3D structure of LOX and its family members. Though few of the previous studies discuss the theoretical structure of LOX the structural orientation of the catalytic site is less discussed and needs to be explored further. In the present study, we have endeavoured to predict a geometrically optimal theoretical structure of LOX and also deduce the proper orientation of the copper ion in the catalytic domain. In the latter part, we have also modelled the copper ion interactions at the copper binding cavity. We have validated the predicted structure and molecular docking studies were conducted to infer the mechanism of its action and binding to the pseudo substrates and experimentally proven inhibitors. It is hoped that the information provided by this study will aid in better understanding on inter and intra molecular interactions of LOX with eventual development of therapeutic applications in LOX mediated diseases.

Methodology: Sequence retrieval and Modelling
In this study, we have employed Ab initio structure prediction approach to model the 3D structure using the ROBETTA server [7], as there are no significant structural templates for homology modelling. Further, this structure was energy minimized and used as a template to generate 1000 models with different conformations using MODELLER 9.10 [8].
Among these 1000 models, the best model with a significant QMEAN score [9], and with no residues in the disallowed regions of Ramachandran plot was selected. Further, this model was fixed for protonation states of histidine residues and its orientation was flipped using Maestro 9.3 (Maestro, Version 9.3, Schrödinger, LLC, New York, NY, 2012) in order to achieve optimal geometry. Similarly, asparagine and glutamine residues were also flipped to achieve optimal geometry.

Copper co-ordination Modelling
The optimal model was visualized in Maestro 9.3 for observing the orientation of histidine residues (292, 294 and 296) spanning the conserved copper binding site; as documented by previous studies [10]. It was observed that the proximity of the imidazole rings favored copper binding. Hence, we implemented constrained molecular dynamics simulation coupled with multiple cycles of energy minimization by OPLS 2005 force field towards achieving the orientation favoring copper ion binding in LOX. Further, we summed up the Cartesian coordinates of copper binding histidine atoms (292N δ , 294N € and 296N € ) and the mean average for X, Y and Z positions were assigned as Cartesian coordinates for the Cu 2+ ion as described in the following equation: (Please see supplementary material for equation 1

and explanation)
The Cartesian coordinates for Y 0 and Z 0 were also similarly derived and assigned to the copper ion. A water molecule was added to the copper ion to satisfy the valence and also to achieve the tetrahedral symmetry as discussed by Ryvkin et al., [11]. Further, the copper ion fixed model was subjected to bond length analysis in order to validate the permitted range of distance (1.9 -2.1 Å) [12].

Model Validation and Refinement
The geometry of the copper ion fixed model was assessed for stereo chemical qualities through PROCHECK [13] and 3D check validation servers [14]. Overall protein topology and domain architecture were also analyzed using the PDBsum server [15].

Molecular Dynamics (MD) Simulation of Homology Models
To infer the stability of the structure predicted, constrained MD simulations were carried out using the Desmond package (Desmond Molecular Dynamics System, version 3.1, D. E. Shaw Research, New York, NY, 2012; Maestro-Desmond Interoperability Tools, version 3.1, Schrödinger, New York, NY, 2012) with an inbuilt OPLS 2005 (Optimized Potentials for Liquid Simulation) force field. As an initial step, the system was prepared for simulation using a predefined water model (simple point charge, SPC) as a solvent in a cubic box with 18 Å x 18 Å x 18 Å dimension as periodic boundary condition. Further, the system was neutralized by adding two Na + counter ions and energy minimized. Finally, the production run was initiated under NPT ensemble conditions for 4 nano seconds. The temperature was set to 300K and maintained throughout by implementing Nose-Hoover thermostat [16] with the pressure set to 1 atm and maintained through Martyna-Tobias Klein pressure bath [17]. Smooth Particle Mesh Ewald method [18] was applied to analyze the electrostatic interactions with a cut-off value of 9.0 Å distance. The Cu 2+ ion, histidine residues involved in Cu 2+ ion interactions and the water molecule bound to Cu 2+ ion were completely constrained during the simulation process. The trajectory sampling was done at an interval of 1.0 pico seconds [16].

Electrostatic potential calculations and Binding pocket prediction
Illustration of the charge distributions of molecular structures is typically rendered through electrostatic potential maps. These maps aid in the identification of sites within the structure to facilitate molecular recognition. The electrostatic interactions between the molecules are generally resolved by the classical Poisson-Boltzmann (PB) equation. In this study, the potential surface for the copper ion fixed model was generated by implementing PB equation through Schrödinger maestro interface. Further, the active site residues were predicted using CASTp server [19]. Finally, contour map for the modelled protein was generated and analyzed for hydrophobic and hydrophilic regions spanning the active site.
The Induced Fit Docking (IFD) of the small molecules on to the predicted active cavity atoms was performed using Schrödinger suite. The final protein model with optimal geometry was imported into Maestro 9.3. Here, the atoms of the active cavity residues (predicted by CASTp) were set as flexible and were assigned as binding site for grid box generation. The ligands were prepared using LigPrep and were docked to the receptor by soften-potential docking with van der Waals radii scaling of 0.70 Å. The resulting 20 best docked conformations with at least one atom within the distance of 5 Å were selected and subjected to geometry optimization, conformational searches and energy minimization. The active cavity residues beyond the range of 5Å, in terms of ligand interactions were set as rigid and those within the 5Å range were set as flexible. Further, the 20 best ligand poses conformations sampled in the initial docking step were re-docked on to the flexible residues within the range of 5Å as followed above. This re-docking was performed using Glide (Extra Precision) XP by soften-potential docking with van der Waals radii scaling of 1.0. Finally, docking score based on OPLS 2005 force field was used to infer the binding affinity of selected small molecules to the receptor. Additionally, Molecular Mechanics Generalized Born Surface Area (MMGBSA) was also calculated to measure the binding free energy (ΔG bind ) of small molecules to the LOX model using Prime/MM-GBSA method [24]. Binding free energy was calculated using the equation: (Please see supplementary material for equation 2 and explanation)

Results: Modelling and Refinement
The initial LOX 3D structure was predicted using ab initio method implementing ROBETTA. Further, a total of one thousand models with varied conformations were generated using MODELLER 9.10 with the initial LOX 3D structure as template. All the models were validated for Ramachandran plot and QMEAN score. Among these, the top ranking model with 92 % of residues in favored region of Ramachandran plot ( Figure 1) and a significant QMEAN score of 0.602 was chosen as the best model.  Then, best model was subjected to refinement by rectifying stereo chemical errors using Schrodinger suite. The refined model was visualized for orientation of histidine residues at the copper binding site and multiple steps of manual minimization by OPLS 2005 were performed to orient these residues such that the copper can form coordinated covalent bonds from His 292N δ , His 294N € and His 296N € favoring tetrahedral symmetry (Figure 2). As a next step, the copper ion was placed in between the interatomic space as per the method discussed by us in the copper co-ordination modelling section. The optimal metal geometry was validated and valence was fixed by adding a water molecule using chimera tool [25]. Finally, the copper fixed LOX model was subjected to multiple steps constrained energy minimization using prime with OPLS2005 as force field. The resultant model was found to have coordinated covalent bond lengths within the allowed distance (1.9 Å-2.2 Å). The copper fixed model was assessed with the Ramachandran plot which showed 90.7 % of residues in favored region with no residues in disallowed region (Figure 3) and overall secondary structure topology analyzed by PDBsum (Figure 4). Further, this model was subjected to constrained molecular dynamics simulation for 4 ns wherein, the copper and its interacting residues were constrained throughout the simulation process. The potential energy of the protein after 4 ns simulation was -307193.771 kJmol -1 . The RMSD trajectory stabilized about 5.5 -6.0 Å after 2 nano seconds of simulation and did not increase significantly after 2 nano seconds. This indicates that the system has evolved into a stable state and has reasonably converged over the production run ( Figure 5A). The radius of gyration analysis showed 0.86 Å of deviation inferring improved relaxation and structural stability of the modelled protein. RMSF graphs also suggest higher flexibility at the Cterminal with few residues showing higher degree of fluctuation ( Figure 5B).

Electrostatic potential graph
Electrostatic potential surface of the copper fixed LOX model was calculated by Poisson-Boltzmann equation and was visualized in Maestro. The charge distribution in the active site cavity was found to be profoundly negative charge which would favor the interactions with positively charged substrates.

Induced fit Docking
The copper fixed model was subjected to induced fit docking with DAP and a group of 3 inhibitors (βAPN, HCys and HCTL) for LOX. The docked complexes were analyzed for docking score, MMGBSA and molecular interaction maps. The substrate and inhibitors were docked to active cavity and the resulting conformations with significant docking score were chosen as probable binding mode (Figure 6). In the docking results, DAP was to found to bind with LOX with a significant docking and MMGBSA score of -7.511kcal/mol and of -41.381 kcal/mol, respectively. DAP also forms hydrogen bonds with Asp 170 , Asp 169 , Asp 353 and Cys 351 of LOX.
In comparison with DAP, the order of binding efficiency for the three inhibitors as follows Hcys > βAPN >HCTL, based on our docking and MMGBSA score Table 1 (see supplementary material). Our docking results indicate that HCys can be an efficient inhibitor for the LOX enzymatic activity. These results also reveal that aspartic acid residues spanning the active cavity region may play the key role in small molecule interactions as it was observed to be a major contributor in hydrogen bonding interactions to all the ligands studied.

Discussion:
In this study, we have attempted to predict the 3D structure of human mature LOX by Ab initio method and also fixed the optimal coordinated covalent interactions of Cu 2+ at the copper binding region. The predicted structure was found to be a valid model as per SAVeS evaluation (http: //nihserver.mbi.ucla.edu /SAVES/). The Root Mean Square Deviation (RMSD) of protein backbone and RMSF (Root Mean Square Fluctuation) of individual residues sampled at periodic intervals during the MD simulation were plotted against the time scale to assess the stability of the model. Here, RMSD plot showed backbone displacement within a range of 0.5 Å after 2 ns and maintained till the end of production run which suggests the stability of the predicted model. Moreover, the radius of gyration plot also infers the compactness of the model. The N -terminal region (169 -220 amino acid residues) of the modelled LOX formed random coil as it was rich in helix breaking residues namely tyrosine and proline. The C -terminal of the modelled LOX was found to form the structural topology with beta strands as similar to that N -terminal of cytokine receptors, which corroborates with the earlier report [1] (Figure 7).
Generally, in copper amine oxidases like LOX, the copper ion plays a crucial role in the structural stability and also in catalytic activity. This region was found to occur in the buried core of most of the crystal structures of copper amine oxidases. Similarly, the copper binding site in the modelled LOX was also found to be in the buried region. Moreover, the copper binding site forms a talon shaped cavity in the LOX model synchronizing to the documented reports [6]. Additionally, we have also modelled the optimal coordinated covalent interactions of Cu 2+ with the histidine residues at the copper binding site with allowed bond lengths ranging from 1.9 -2.2 Å and a water molecule interaction in the tetrahedral geometry which corroborates with the reported crystal structure of E.Coli amine oxidase [12].
In amine oxidases, the substrate binding and positioning are mainly guided by charge -charge interactions. These substrates are generally cationic and found to interact with the key anionic residues at active cavity. Similarly, the active cavity of the modelled LOX structure was observed to be anionic. Further, the proton abstraction is usually catalyzed by aspartic acid in case of all amine oxidases [26,27]. In this study, the IFD results for LOX with pseudo substrate DAP and all the inhibitors also strongly infer Asp residues as major contributors of hydrogen bonding interactions. All these findings strongly suggest the plausibility and reliability of the modelled structure and it inter molecular interactions. Till date, the paucity of information on structural aspects of LOX remains as a limiting factor for understanding its role in cellular processes. Hence, the outcomes of this study shall evoke new dimensions towards exploring structure-function relationships of LOX.

Conclusion:
In this study, we predicted the optimal structure of LOX with coordinated covalent orientation of copper ion at its catalytic cavity. IFD was performed to understand the molecular interactions of LOX with its modulators which inferred aspartic acid residues as the key contributors towards intermolecular interactions. The predicted structure was validated by stereo chemical checks, MD studies and leads from the literature. All these findings reinforce the higher plausibility of the predicted structure and its intermolecular interactions. The insight of this study will pave way for design and development of novel therapeutic molecules potentially modulating LOX activity.