Identification and Evaluation of Inhibitors of Lipase from Malassezia restricta using Virtual High-Throughput Screening and Molecular Dynamics Studies

Recent studies revealed the role of lipase in the pathogenicity of Malassezia restricta in dandruff and seborrheic dermatitis (D/SD). The lipase from M. restricta (Mrlip1) is considered a potential target for dandruff therapy. In this work, we performed structure-based virtual screening in Zinc database to find the natural bioactive inhibitors of Mrlip1. We identified three compounds bearing superior affinity and specificity from the Traditional Chinese Medicine database (~60,000 compounds), and their binding patterns with Mrlip1 were analyzed in detail. Additionally, we performed three sets of 100 ns MD simulations of each complex in order to understand the interaction mechanism of Mrlip1 with known inhibitor RHC80267 and the newly identified compounds such as ZINC85530919, ZINC95914464 and ZINC85530320, respectively. These compounds bind to the active site cavity and cause conformational changes in Mrlip1. The Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) studies suggested that the average binding energy was stronger in the case of Mrlip1-ZINC85530919 and Mrlip1-ZINC95914464. The selected natural inhibitors might act as promising lead drugs against Mrlip1. Further, the present study will contribute to various steps involved in developing and creating potent drugs for several skin diseases including dandruff.


Introduction
Lipases (E.C. 3.1.1.3) belong to the hydrolase family that act on carboxylic ester bonds [1]. Monoacylglycerol (MAG) and diacylglycerol (DAG) lipases have been taken in consideration because of their physiological functions as a biocatalyst [2] and applications in oil and fat modification [3]. Lipases are one of the most widely used classes of enzymes in biotechnology and synthetic organic chemistry [4]. In spite of the fact that Malassezia restricta lipase (Mrlip1) is a mono-and diacylglycerol lipase secreted from M. restricta KCTC27527, it was isolated and identified from a patient with dandruff and seborrheic dermatitis (D/SD) [5,6]. Dandruff is a common scalp disorder that has a prevalence of nearly 50% in the worldwide population [7]. Several studies have claimed associations between M. restricta and specific skin diseases [8]. The Malassezia spp. lipases degrade the triglycerides of human sebum and consume specific saturated fatty acids which in turn cause irritation in individuals' skin with dandruff and seborrhoeic dermatitis [9].
It has been well established that pathogenic fungi produce extracellular lipases to breach the host tissue barrier and enable them to penetrate the tissue. The examples include lipases of Candida albicans that contribute to yeast-hypha transition, colonization and invasion into the host tissues [10,11]. Moreover, the lipases were also identified in M. furfur and M. globosa, and their biochemical characteristics were analyzed [9,12].
Current treatments of D/SD are mostly anti-fungal agents such as zinc pyrithione, ketoconazole, coal tar, selenium sulfide, and general lipase inhibitors [13]. However, efficient, stable and non-polymeric therapeutic agents are much needed. The activity of Mrlip1 is likely to be maximal under the conditions that M. restricta typically faces at the host skin surface and involved in dandruff diseases [6].
The present work aims to seek a new class of inhibitors specifically targeting the Mrlip1 by Virtual High-throughput Screening (vHTS) and molecular dynamics (MD) simulation methods. We have performed structure-based vHTS of about 60,000 small compounds in Traditional Chinese Medicines (TCM)-based naturally occurring compounds from non-commercial ZINC database followed by several sets of 100 ns all-atom MD simulations to find out potent inhibitors of Mrlip1. The results suggested that ZINC95914464 is a potent bioactive compound that represents novel hits that could serve as the starting point for the development of more potent anti-dandruff therapeutic agents.

Results and Discussion
The Malassezia restricta lipase (Mrlip1) is mainly involved in dandruff progression [6]. Amongst all species of dandruff causing Malassezia, the most dominated species is M. restricta [7]. Therefore, a total of about 60,000 small compounds were screened against the structure of Mrlip1, and 80 top-ranked compounds that have the highest binding affinity were selected for further screening. All novel hits were accurately fitted within the active site of Mrlip1 and were further evaluated for drug-likeness using various tools.

vHTS: Molecular Docking
The screening of compounds library produced log-files and output-files, which contain binding affinity scores and docked poses for individual compounds in the library. These log-files and output-files were subjected to a screen-out based on their binding affinities, docking score and binding orientation for Mrlip1. The number of natural compounds having a good binding affinity score were selected further in the search for potential inhibitors of Mrlip1.

Hit Selection and Drug-Ability Assessment
Initially, the compounds were filtered out to get the highest binding affinity natural compounds from the 60,000 screened-compounds, which were extracted by a python script. We obtained the 80 highest binding affinity natural compounds (Table S1). Further, these compounds were subjected to further screening based on their physicochemical properties, where 25 compounds were qualified in specific cut-off values of drug-likeliness. The compounds were selected based on parameters such as hydrogen bond donors less than 5, hydrogen bond accepters less than 10, rotatable bonds less than 10, molecular weight less than 750 Dalton, and logP less than 10 (Table 1).
However, some compounds violated the Lipinski rule of five like many FDA approved drugs do, as they have a molecular weight of >500 Dalton and logP value >5, but this is acceptable [14]. ADMET properties of these compounds were predicted, where 14 compounds were showing acceptable ADMET properties ( Table 2). All these results indicate that the selected compounds show an ideal characteristic behaviour of a drug-like molecule. Furthermore, interaction analysis was carried out to get selective compounds specific to the binding site of Mrlip1, and finally, the three best compounds were selected and one additional reported compound was taken as a control inhibitor RHC80267 (Table 3). These selected compounds passed the Pan-assay interference compounds (PAINS) filters, where no alert of PAINS pattern was found for these compounds, which supports the specificity of the selected compounds towards Mrlip1. Based on our results, we proposed that these three drug-like compounds may be considered as potent inhibitors of Mrlip1, showing appreciable binding affinity and specificity to the Mrlip1 binding pocket, which can reduce the accessibility of the substrate and thus, interfere with enzyme activity.

Structural Analysis of Mrlip1 Complexes
The structural analysis of Mrlip1 complexes suggested that the catalytic pocket consists of several important residues such as Ser 171 , Tyr 54 , Thr 101 , Ile 106 , Thr 107 , Phe 278 , His 281 , Gln 282 , Gln 289 , Ala 292 and Phe 294, which are responsible for its catalytic activity and substrate binding. These residues fall on lid and flap domains and also consider the catalytic triad (Ser 171 , Asp 228 and His 281 ) interactions. The lid protects the active site and hence, is responsible for catalytic activity [15]. It has been predicted that the interaction with these catalytic pocket residues will significantly reduce the catalytic activity of Mrlip1. The selected compounds show close interactions with active site residue Ser 171 and His 281 of Mrlip1 (Figures 1-3).  All these compounds are present in the deep cavity of Mrlip1 and mimick the binding pose with one another (Figure 1). The residues of Mrlip1 such as Ser 171 , Tyr 279 , Gln 282 , Arg 83 , and Thr 101 form hydrogen bonds with the compounds in addition to several van der Waals and other interactions. Apart from these catalytic residues, many important interactions are offered by the catalytic pocket residues such as Tyr 54 , Asp 56 , Asn 77 , Thr 107 , Leu 172 , Phe 278 , His 281 , Gln 289 , Ala 292 , and Phe 294 to the compounds to properly hold in the cavity of Mrlip1 (Figures 2 and 3). Surface representations clearly indicated that the compounds reside in the internal part of the catalytic pocket and specifically bind with Mrlip1 important pocket residues with significant affinity ( Figure 1B).

MD Simulation Data Analysis
To ensure the MD simulation data reproducibility and test the convergence of the result, the final production phase of each system was carried out with three independent MD set runs [16,17]. The replica has been shown as Run 1, Run 2, and Run 3 in Figures S1-S10. Furthermore, a detailed analysis was carried out using all trajectories.

Average Potential Energy of Complex Systems
To ascertain the equilibration of the systems prior to MD analysis, the average potential energy of free Mrlip1, Mrlip1-RHC80267, Mrlip1-ZINC85530919, Mrlip1-ZINC95914464, and Mrlip1-ZINC85530320 were monitored. The constant temperature fluctuations at 300 K for each system suggest the stable and accurate nature of the MD simulations performed. The average potential energy for free Mrlip1, Mrlip1-RHC80267, Mrlip1-ZINC85530919, Mrlip1-ZINC95914464 and Mrlip1-ZINC85530320 were found to be −642,080 kJ/mol, −630,167 kJ/mol, −629,866 kJ/mol, −629,362 kJ/mol, and −632,500 kJ/mol, respectively.

Conformational of Mrlip1
Binding of a compound in the catalytic pocket of the protein can lead to large conformational variations in a protein [18][19][20]. Root mean square deviation (RMSD) is one of the most important fundamental properties for establishing whether the protein is stable and close to the experimental structure [21]. The average RMSD values for free Mrlip1, Mrlip1-RHC80267, Mrlip1-ZINC85530919, Mrlip1-ZINC95914464 and Mrlip1-ZINC85530320 were found to be 0.25 nm, 0.21 nm, 0.27 nm, 0.20 nm, and 0.28 nm, respectively (Table 4). The RMSD plot suggests that the binding of RHC80267 and ZINC95914464 stabilized the Mrlip1 and leads to less structural deviations from its native conformation ( Figure 4A). In the case of the Mrlip1-ZINC85530919 and Mrlip1-ZINC85530320 complex, the binding of ZINC85530919 and ZINC85530320 to the Mrlip1 active pocket showed to have initial low fluctuations until 10-25 ns of the MD trajectories, thereafter it attained a high RMSD value. The orientation of ZINC85530919 and ZINC85530320 in the active pocket of Mrlip1 showed to have the least RMSD values and was found to be equilibrated during the 100 ns MD simulations. The RHC80267 and ZINC95914464 compounds showed continuous fluctuations in the active pocket of Mrlip1, which is possibly due to different orientations ( Figure 4B). Further analyses of each complex was carried out at different time intervals such as 1 ns, 40 ns, 80 ns and 100 ns. The orientations of RHC80267 in the catalytic pocket were totally different mainly due to the molecular structure of the compound. The known RHC80267 inhibitor contains two rings on both sides connected with flexible chains. Therefore, we observed one ring in the catalytic pocket and the other ring on the surface of the catalytic pocket ( Figure 5A-D). On the other hand, ZINC95914464 compound also has a similar molecular structure, but it is not more flexible than RHC80267. These properties of the ZINC95914464 compound may make it a better drug-like molecule, since it occupied the catalytic pocket completely ( Figure 5E-H). The ZINC85530919 and ZINC85530320 compounds showed similar orientation within the catalytic pocket. The catalytic pocket of Mrlip1 is highly selective in nature for the substrates. Therefore, very few interactions were formed by ZINC85530919 and ZINC85530320 compounds in the catalytic pocket ( Figure 5I-L). Due to this reason, they cannot occupy the whole catalytic pocket. The replicas of all complexes showed a similar pattern of binding interaction with Mrlip1 ( Figure S2).  Vibrations around the equilibrium are not random but depend on the local structure flexibility. To calculate the average fluctuation of all residues during the simulation, the root-mean square fluctuation (RMSF) of the Mrlip1 upon compounds binding was plotted as a function of residue number ( Figure 4C). The RMSF plot showed that residual fluctuations are present in Mrlip1 in several regions of the protein structure. These residual fluctuations were found to be minimized upon binding of RHC80267 and ZINC95914464 and maximized with ZINC85530919 and ZINC85530320 at the region spanning the lid and flap domain.
Radius of gyration (R g ) is a parameter linked to the tertiary structural volume of a protein and has been applied to obtain insight into the stability of a protein in a biological system. A protein is supposed to have a higher radius of gyration due to less tight packing. The average R g values for free Mrlip1, Mrlip1-RHC80267, Mrlip1-ZINC85530919, Mrlip1-ZINC95914464 and Mrlip1-ZINC85530320 were found to be 1.67 nm, 1.64 nm, 1.65 nm, 1.64 nm, and 1.67 nm, respectively. R g plot suggested that the Mrlip1 attained more tight packing in Mrlip1-RHC80267 and Mrlip1-ZINC95914464, and when bound to Native, ZINC85530919 and then ZINC85530320 ( Figure 4D).

Hydrogen Bonds Analysis
Hydrogen bonding between a protein and compounds provides a directionality and specificity of interaction that is an important aspect for molecular recognition [23]. In order to validate the stability of the docked complexes, the hydrogen bonds paired within 0.35 nm between Mrlip1 and compounds in Mrlip1-RHC80267, Mrlip1-ZINC85530919, Mrlip1-ZINC95914464 and Mrlip1-ZINC85530320 were calculated in a solvent environment during the MD simulations. We found that RHC80267 and ZINC95914464 strongly bound to the active pocket of Mrlip1 with 2-6 hydrogen bonds, while ZINC85530919 and ZINC85530320 bound to the active pocket of Mrlip1 with 1-3 hydrogen bonds with least fluctuations (Figure 7).

Secondary Structure Changes upon Ligands Binding
The purpose of this analysis is to measure the secondary structure content of a protein as a function of time. The secondary structure assignments in protein such as the α-helix, β-sheet and turn were fragmented into individual residues for each time step. The average number of residues that participated in secondary structure formation in the case of complexes was slightly decreased due to an increase in the percentage of coils and a decrease in the β-sheet in comparison with Mrlip1 (Table 5). There was a slight decrease in the percentage of the β-sheet, and an α-helix was found in the case of Mrlip1-ZINC85530919 (Figure 8).

Principal Component Analysis
The principal component analysis (PCA) or essential dynamics (ED) reflects the overall expansion of a protein during different simulation [24]. In this method, the dynamics of Mrlip1 was calculated using gmx covar module with respect to the backbone. PCA identifies the large scale average motion of a protein, thus revealing the structures underlying the atomic fluctuations [25]. The sum of the eigenvalues is a measure of the total motility in the system. It can be used to compare the flexibility of a protein under different conditions [26]. The trace of the covariance matrix and eigenvalues for free Mrlip1, Mrlip1-RHC80267, Mrlip1-ZINC85530919, Mrlip1-ZINC95914464, and Mrlip1-ZINC85530320 were found to be 131.214 nm 2 , 146.472 nm 2 , 151.465 nm 2 , 212.819 nm 2 , and 196.618 nm 2 , respectively. The eigenvalues for free Mrlip1 were found to be low in comparison with complexes that clearly indicate that a random fluctuation increases upon compounds binding; and Mrlip1-ZINC95914464 showed higher eigenvalues, which indicates a higher expansion of Mrlip1 i.e., low compactness of complexes at 300 K. The multidimensional atomic covariance matrix for each atom pair covariance is depicted (Figure 9). The gmx anaeig module reads a set of eigenvectors and eigenvalues as input files and returns to project an MD trajectory along selected eigenvector values. The 2D projections of trajectories on eigenvectors showed overlap between free Mrlip1, Mrlip1-RHC80267, and Mrlip1-ZINC95914464 complexes. The results also suggested that the difference in the position of atoms is due to the binding of complexes with Mrlip1. The root mean square atomic fluctuations of Mrlip1 were also recorded during the PCA calculations. The Eigenvector components were further resolved into x, y and z directions.

Gibbs Free Energy Landscape
The Gibbs free energy landscape was also calculated using gmx covar, gmx anaeig, and gmx sham using the projections of their own first (PC1) and second (PC2) eigenvectors, respectively. The color-coded Gibbs free energy landscape is also depicted in Figure 10. Gibbs free energy landscape inspects the direction of the fluctuation in the two systems for all Cα atoms of the free Mrlip1, Mrlip1-RHC80267, Mrlip1-ZINC85530919, Mrlip1-ZINC95914464, and Mrlip1-ZINC85530320 complex structures from trajectories. The corresponding free energy contour map with a deeper blue color indicates a lower energy. It has been found that the main free energy well in the global free energy minimum region was completely changed after binding of these compounds. These free energies well indicate the stable conformational states of this molecule. The binding of these compounds to Mrlip1 leads to different global minima of Mrlip1 during the 100 ns MD simulations. The energy landscape exhibits several clearly distinguishable minima, which correspond to the metastable conformational states that are separated with a small energy barrier. The most metastable conformational states were seen in the binding of Mrlip1-RHC80267, Mrlip1-ZINC85530919, and Mrlip1-ZINC95914464, in which local minima distributed to about three to four regions within the energy landscape. The Free Mrlip1 and Mrlip1-ZINC85530320 formed just two to three metastable conformations during the whole population of trajectories ( Figure S10).

MMPBSA Binding Energy Analysis
The extraction of binding energies from the stable short region of 40 to 50 ns simulation, 200 trajectory frames were used for the mmpbsa calculations using polar and apolar solvation parameters. The aim of this analysis is to establish the energies associated with the binding of compounds to Mrlip1 during the MD simulations. The average binding energy for Mrlip1-RHC80267, Mrlip1-ZINC85530919, Mrlip1-ZINC95914464 and Mrlip1-ZINC85530320 were found to be −98.51 kJ/mol, −232.28 kJ/mol, −183.17 kJ/mol, and −85.71 kJ/mol, respectively (Table 6).
(visual molecular dynamics) [29] and Discovery studio visualizer [30] were used for visualization purposes. In addition, the online web-servers and data repository resources like NCBI, PDB, ZINC database, PubChem, SwissADME [31], and CarcinoPred-EL [32] were used for retrieval, evaluation and analysis of the data.

Preparation of Target Protein and Natural Compounds Library
The three-dimensional (3D) structure of open conformation of Mrlip1 was obtained from our previous work [33]. The open conformation of Mrlip1 was validated by well-studied Rhizomucor miehei lipase (PDB ID: 4TGL) which was co-crystalized with DEP (Di-Ethyl-Phosphonate) in its catalytic pocket. The same pose and orientation of DEP was modelled in Mrlip1 lipase and used for vHTS. A library of natural compounds based on the TCM database containing 60,000 small compounds was obtained from non-commercial ZINC database [34]. The library contains processed chemical structures of TCM natural compounds in 3D file formats which were used in this study to perform structure-based vHTS.

Hit Selection and Drug-Ability Evaluation
Based on binding affinities and scoring, the top 80 hits were selected, which showed higher binding-affinity towards Mrlip1. The selected compounds were further analyzed based on physicochemical and ADME properties to find out safe and effective drug-like compounds. These properties including toxicity and carcinogenicity were predicted using web-based applications such as Swiss-ADME [31], PreADMET and CarcinoPred-EL [32]. Further, the selected compounds were checked for a PAINS (Pan-assay interference compounds) pattern using the SwissADME web server, accessed on: 23 June 2018) [35] and ZINC 15 chemistry pattern database, accessed on: 25 June 2018) [36] and then further interaction examinations were carried out to avoid a false positive and get selective compounds with high specificity towards the binding pocket of Mrlip1.

Mrlip1 Structural Analysis: Visualization and Evaluation
The visual inspection of docked conformations of compounds with Mrlip1 was performed using PyMOL and Discovery Studio visualizer LigPlot + visualization tools. These softwares create high quality animated 3D as well as 2D figures of Mrlip1 and chemical compounds. Apart from the visual inspections, numerous measurements such as bond length, distance between residues, and distance between the Mrlip1 and compounds were calculated.

MD Simulations
Molecular dynamics (MD) simulation is a convenient method to understand the chemistry of biological macromolecules at the molecular level [37,38]. Many mysterious biological functions in proteins and their profound dynamic mechanisms can be revealed by studying their internal motions [39][40][41][42]. MD simulations were performed on free Mrlip1, Mrlip1-RHC80267, Mrlip1-ZINC85530919, Mrlip1-ZINC95914464, and Mrlip1-ZINC85530320 at 300 K at the molecular mechanics level implemented in the GROMACS 5.1.2 [43] using the GROMOS96 43a1 force-field [44]. The compounds topology and force-field parameters were generated using the PRODRG server [45]. Charges in the topology file were manually corrected. The additional compound atoms were merged in the complex topology files, and the parameters for all the compounds were included in the system topology.
The free Mrlip1, Mrlip1-RHC80267, Mrlip1-ZINC85530919, Mrlip1-ZINC95914464, and Mrlip1-ZINC85530320 were soaked in a cubic box of water molecules with an initial dimension of 8 nm, i.e., setting the box edge 1 nm from the molecule periphery using the gmx editconf module for creating boundary conditions and gmx solvate module for solvation. The Simple Point Charge (spc216) water model was used to solvate the protein. The charges on the free Mrlip1 and its complexes were neutralized by the addition of Na + and Cl − ions using the gmx genion module to maintain neutrality, preserving a physiological concentration (0.15 M). Energygrps in the molecular dynamics parameters (mdp) file was used to examine the interactions between Mrlip1 with compounds and the particle-mesh Ewald method [46] was applied. Then, the MD system was minimized by means of steepest descent (1500 steps). The temperature was then raised (0 to 300 K) during the equilibration period of 100 ps under periodic boundary conditions at a constant volume.
Equilibration was completed in two constant stages such as NVT (number of particles, volume, and temperature) ensemble and NPT (number of particles, pressure, and temperature at 100 ps) ensemble. The final production phase (100 ns) was achieved thereafter at 300 K. Several Gromacs modules were used to analyze the MD trajectories such as gmx energy, gmx rms, gmx confirms, gmx rmsf, gmx gyrate, make_ndx, gmx hbond, gmx anaeig, gmx sham, gmx do_dssp, gmx_covar, and gmx sasa. The graphical presentation of the 3D models was prepared using VMD (Visual Molecular Dynamics) [29] and PyMOL.

MMPBSA Calculation
Molecular mechanics Poisson Boltzmann surface area (MMPBSA) is an attractive approach which has been used effectively to reproduce and justify experimental findings and to improve the results of virtual screening and docking [47]. A short MD trajectory was extracted from the stable region of each complex for MMPBSA estimations [48]. The binding energy was calculated using MMPBSA approaches of the g_mmpbsa package [49]. It was calculated using the following equation: where, G Complex signifies the total free energy of the binding complex, and G Protein and G Ligand are the measures of the total free energies of the individual protein and compounds, respectively.

Conclusions
Identifying potential and specific inhibitors for Mrlip1 is a promising approach to address skin diseases. In this study, the Traditional Chinese Medicine database has been screened against Malassezia restricta lipase (Mrlip1) to find promising and highly effective anti-dandruff inhibitors. The possible inhibitors of Mrlip1 such as ZINC85530919, ZINC95914464 and ZINC85530320 were selected on the basis of high binding affinities, interactions with catalytic triad residues, and predicted activity scoring values. Based on several screenings, three best hits were identified and are expected to be possible inhibitors of Mrlip1. Various examinations such as drug-likeliness, ADME and toxicity were applied on all three compounds. All selected compounds interact with catalytic pocket residues and occupy the same binding pocket as RHC80267, a well-known diacylglycerol lipase inhibitor. These compounds have a better average binding affinity. The MD simulation study revealed that these compounds formed stable conformations within the catalytic pocket and exhibit different conformational changes in the Mrlip1. Amongst all the identified inhibitors, ZINC95914464 was found the best stability in the catalytic pocket. Therefore, ZINC95914464 can be a new inhibitor that could serve as the starting point for the development of more promising anti-dandruff therapeutic agents.