Computational study of potential inhibitors for fat mass and obesity-associated protein from seaweed and plant compounds

Protein and ligands preparation

FTO protein with 3LFM as corresponding PDB ID was retrieved for this study (Han et al., 2010). It possessed chain A that corresponds to the FTO protein and consisted 495 amino acids from Homo sapiens. Some missing amino acids were found to be in the loops which showed no interaction at the active site and were excluded from the study. The entire structure had a resolution of 2.50 Å. Prior to docking studies, ligand molecules in the protein structure, non-interacting ions and water molecules were removed using BIOVIA Discovery studio client 2020 (Dassault Systèmes, BIOVIA, Discovery Studio Modelling Environment, Release 2017, San Diego: Dassault Systèmes, 2016).

A library of 100 plant compounds having anti-obesity activity for the FTO protein was retrieved from published literature. The plant sources of these compounds are mentioned in Table S1 and their corresponding structures were downloaded from PubChem database. In addition, 1,110 seaweed compounds retrieved from Seaweed Metabolite Database were provided in Table S2 (Davis & Vasanthi, 2011). The 3D structure of all the ligand molecules were generated using Corina and converted to PDB format using an Open Babel module (O’Boyle et al., 2011).

Virtual screening

PyRx v0.8 (Trott & Olson, 2009) is a structure-based virtual screening tool with an inbuilt module of Open Babel for ligand processing and AutoDock Vina 1.12 (Dallakyan & Olson, 2015) for molecular docking. The 3D structures of the ligands and the protein were loaded in PyRx. The conjugate gradient algorithm was implied to these structures for energy minimization using Universal Force Field (UFF). Energy minimized structures were then converted to PDBQT format. The docking site was defined by 3-methylthymidine (DT) from the same 3LFM structure data. Exhaustiveness modes were set to eight and a total of eight different docking poses were generated and analysed based on binding affinity and interactions with the protein.

Drug-likeness and oral toxicity prediction

The selected molecules were subjected to ADME predictions using SwissADME (Daina, Michielin & Zoete, 2017) to analyse drug likeness property. In an attempt to find optimal bioavailability of the drug, Lipinski’s rule of five (RO5) (Lipinski et al., 2001) was used as a filter for the virtual screening of compounds. Furthermore, with the aim of achieving significant oral non-toxic drugs, toxicity of the selected molecules was evaluated using ProTox-II server (Banerjee et al., 2018). ProTox-II predicts the median lethal dose (LD50) values in mg/kg body weight (Drwal et al., 2014) based on the Globally Harmonized System of classification of labelling of chemicals (Winder, Azzi & Wagner, 2005).

Molecular docking

The in silico analysis for the selected compounds with 3LFM were performed using AutoDock 4.2.6 with ADT tools of 1.5.6 version (Morris et al., 1998). The target protein 3LFM was refined by removal of hetero atoms and water molecules followed by addition of polar hydrogen. Kollaman charges were computed for the receptor and saved in PDBQT format. The calculation of Gasteiger charges was subjected to energy minimized conformations of the ligands. The torsion counts of rotatable bonds were further defined for the ligands and saved in PDBQT format. The Lamarckian Genetic algorithm parameters were applied with a population size of 300, mutation rate of 0.02, energy evaluations of 2.5 × 10⁶, crossover rate of 0.08 and 100 long GA runs. The docked conformations of each protein and ligand were predicted with the energy values in kcal/mol. The docked output of the protein-ligand complex was rendered using UCSF Chimera v1.15 software (Pettersen et al., 2004). Based on the binding affinity, the best docked complexes were visualized using BIOVIA Discovery studio client 2020 (Dassault Systèmes, BIOVIA, Discovery Studio Modelling Environment, Release 2017, San Diego: Dassault Systèmes, 2016). The efficiency of the ligand molecules was assessed using the ligand efficiency calculation as stated by Hopkins, Groom & Alex (2004) LE represents average binding energy in relation to number of heavy atoms in a molecule.

Molecular dynamics (MD) simulation study

The molecular dynamics simulation study was performed using GROMACS 2020.1 (GROningen Machine for Chemical Simulations) software (Hess et al., 2008) to explore the structural stability of protein and protein-ligand complexes. The ligand topology and the parameter files were generated using PRODRG server (Schüttelkopf & van Aalten, 2004) and protein topology was generated using GROMOS96 53a6 force field (Oostenbrink et al., 2004) to produce energy minimized conformation of all protein-ligand complexes. These complexes were solvated using the Extended Simple point charge (SPC-E) water cubic box model (Wu, Tepper & Voth, 2006). Protonation state of the residues were determined using PROPKA v2.0 under physiological pH 7.4 (Bas, Rogers & Jensen, 2008) and 14 Na⁺ ions was added to each system to maintain the overall neutrality. Further energy minimization was performed using 50,000 steepest descent steps for each protein-ligand complex. After energy minimization process, each system was equilibrated under NVT (number of particles, volume and temperature) followed by NPT (number of particles, pressure and temperature) for 100 ps time scale of position restrain. Constant temperature was maintained at 300 K using V-rescale temperature coupling method and pressure was maintained at 1 atm with Parrinello-Rahman coupling method (Martoňák, Laio & Parrinello, 2003). The particle-mesh ewald (PME) (Toukmaji et al., 2000) and Verlet cutoff scheme at 10 kJ mol⁻¹ (Grubmüller et al., 1991) were used to measure long range electrostatic interactions and van der Waals interactions, respectively. The lincs (Linear Constraint Solver) algorithm was applied for covalent bond constraints (Hess et al., 1997). Finally, MD simulation step was carried out to analyse the stability of protein-ligand complexes for 100 ns time scale. The trajectories generated by the simulation were analysed using in-built GROMACS script to evaluate RMSD (root mean square deviation), radius of gyration (Rg), RMSF (root mean square fluctuation) and hydrogen bonds (H-bond) formation. The graphs were generated and visualized using XMGrace tool.

Binding free energy calculation

Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) was used to calculate the binding free energy (ΔG) of protein-ligand complexes using the g_mmpbsav5.12 package in the GROMACS platform (Kumari, Kumar & Lynn, 2014). Furthermore, to understand the contribution of each amino acid and its dynamic behaviour, the binding free energy which includes the potential energy (electrostatic and van der Waals’ interactions) and free solvation energy (polar and nonpolar solvation energies) of each protein-ligand complex were estimated using the MM/PBSA method. Finally, the stable last 20 ns trajectories of MD simulation for each protein-ligand complexes were selected to compute binding free energy.

Virtual screening

The virtual screening of 1,210 compounds was performed using AutoDock Vina in PyRx. Based on the binding energy and their interactions, 443 compounds with binding energy lower than −7.00 kcal/mol were selected for further analysis. The binding energies of 443 compounds are shown in Table S3. The docking protocol was validated by docking 3-methylthymidine ligand into the active site of 3LFM.

Drug-likeness

Drug-likeness of the screened compounds were evaluated using RO5 parameters. Finally, out of 443 compounds, 369 compounds satisfied Lipinski’s properties and exhibited drug-likeness indicating possible permeability and absorption of molecules through biological membranes. RO5 parameters for the 369 compounds are provided in Table S4.

Toxicity prediction

Acute oral toxicity was evaluated based on ProTox II results. According to Winder, Azzi & Wagner (2005), 30 compounds exhibited non-toxicity with an LD50 value greater than 5,000 mg/kg under class VI. These non-toxic compounds were chosen for the next step of analysis. The toxicity class with LD50 values in mg/kg is listed in Table S5.

Docking studies

Molecular docking studies were carried out between FTO protein and selected 30 ligands using AutoDock 4.2.6 tools to predict best active compounds. The results obtained proposed that out of 30 compounds docked with the FTO protein, four marine-derived compounds BT012 (24ε-hydroperoxy-6β-hydroxy-24-ethylcholesta-4,-28(29)-dien-3-one), BD064 ((6β,24ε)-24-Hydroperoxy-6-hydroxystigmasta-4,28-dien-3-one), RL442 (10-acetoxyangasiol) and RL074 (2-hydroxyluzofuranone B) exhibited high binding affinity values with more negative docking scores when compared to the reference molecule 3-methylthymidine (DT) as well as plant derived compounds (Table 1). The affinity of the ligands in the studied was assessed by ligand efficiency (Table 1). All intermolecular energy interactions and key residues involved in the interactions of 30 compounds are summarized in Table S6.

Table 1:

Binding energy and interactions of ligands docked with FTO protein (3LFM).

Compound ID	Compound name	Chemical structure	Binding energy (kcal/mol)	Ligand efficiency (kcal/mol)	Interactions with amino acid residues of 3LFM
					Hydrogen bond interactions	Hydrophobic interaction	Pi interactions
3-methylthymidine (DT)	FTO inhibitor		−6.61	0.37	Arg96, Tyr106, Asn205, Arg322	Thr92, Pro93, Val94, Tyr108, Leu203, Ser229, Asp234, Glu234, Thr320	Val228, His231
BT012	24ε-hydroperoxy-6β-hydroxy-24-ethylcholesta-4,-28(29)-dien-3-one		−9.12	0.30	Thr92, Tyr106, Glu234	Pro93, Tyr108, Leu109, Val228, His231	–
BD064	6β,24ε)-24-Hydroperoxy-6-hydroxystigmasta-4,28-dien-3-one		−8.41	0.36	Pro93, Tyr106, His231, Glu234	Val94, Val228	His231
RL442	10-acetoxyangasiol		−8.01	0.31	Pro93, Lys216, Ser229	Tyr108, Leu109	–
RL074	2-hydroxyluzofuranone B		−7.5	0.35	Pro93, Asp233, Glu234	Val94, Val228	His231

DOI: 10.7717/peerj.14256/table-1

The docking results obtained pointed out that the reference molecule DT had four conventional hydrogen bonds with Arg96, Tyr106, Asn205 and Arg322 have a binding energy of −6.61 kcal/mol. Nine residues are involved in the hydrophobic interaction and Val228 and His231 residues are involved in Pi alkyl formation (Fig. 1A).

BT012 have binding energy value of −9.12 kcal/mol. Two conventional hydrogen bonds with Tyr106 and Glu234 and a single carbon-hydrogen bond with Thr92 were observed in the binding of BT012 to the active site of 3LFM. In addition, alkyl interaction with Leu109, Val228 residues and Pi alkyl interactions with Pro93, Thr108 and His231 were also observed and are represented in Fig.1B. The conformational binding of BD064 exhibited binding energy of −8.14 kcal/mol which possessed three conventional hydrogen bonds with Tyr106, His 231, Glu 234 and a single carbon-hydrogen with Pro93. Alkyl interactions of Val94 and Val228 were observed in addition to the Pi lone pair of His231 residue (Fig. 1C).

The docked results from compound RL442 had two conventional hydrogen bonds with Lys216 and Ser229 and one carbon-hydrogen bond with Pro93. Hydrophobic interactions formed with Tyr108 and Leu109 had a binding affinity of −8.01 kcal/mol (Fig. 1D). The interactions of the compound RL074 with the active site of 3LFM is depicted in Fig. 1E. The conformational binding exhibited binding energy of −7.5 kcal/mol. The obtained docking results pointed out that formation of conventional hydrogen bonds with Asp233, Glu234 and carbon-hydrogen bonds with Pro93 were important residues to maintain energy conformations. Alkyl interactions of Val94 and Val228 residues were formed, in addition to this Pi lone pair of His231 residue was noted (Fig. 1E).

The four selected marine derived compounds such as BT012 from Turbinaria conoides (Sheu et al., 1999), BD064 from Dictyota bartayresiana (Rao et al., 1994), RL442 from Laurencia sp (Vairappan et al., 2010) and RL074 from Laurencia saitoi (Sun et al., 2005) showed a lower docking energy score in comparison with DT as well as with plant derived compounds. Pro93, Val94, Tyr106, Tyr108, Leu109, Glu234, Val228 and His231 residues are found to be in common with the reference molecule. These compounds also exhibited strong interaction with the target protein and may be considered as potential inhibitors of FTO. The docked protein-ligand complexes were eventually used to study the structural, dynamic as well as binding behaviours of ligand interactions with the active site of 3LFM using MD simulation studies.

Molecular dynamics (MD) simulation

MD simulation trajectories were accomplished to assess the detailed structural and dynamic interaction of the FTO protein along with the FTO-ligand complexes (BT012, BD064, RL442, RL074 and DT) for 100 ns time scale. The RMSD-protein backbone, RMSF analysis, hydrogen bond analysis and Rg scores of complexes during the entire simulation period were analysed.

The RMSD of the backbone atoms were assessed to understand the complex stability of complexes over 100 ns time scale. The RMSD of the trajectories revealed that apo FTO reached stability after 65 ns whereas all the other FTO-ligand complexes took 20–25 ns to reach mere stability throughout the system. The average RMSD values of FTO protein, FTO-BT012 complex, FTO-BD064 complex, FTO-RL442 complex, FTO-RL074 complex and FTO-DT were found to be 0.39, 0.48, 0.41, 0.34, 0.34 and 0.39 nm, respectively. As it was seen that RMSD deviation of protein backbone in selected complexes ranges between 0.2–0.5 nm and does not exceed 0.5 nm which indicates the system to be merely similar and aids in maintaining structural integrity as depicted in Fig. 2. The ligand RMSD plot exhibited stable dynamic profiles in the range of 0.05 to 0.25 nm. The ligand RMSD plot for the ligand complexes are given in Fig. 3. RL074 and RL442 ligands showed stable deviation when compared to reference ligand over 100ns time scale.

The dynamic flexibility and mobility of the residues were computed using RMSF calculations for FTO protein and FTO-ligand complexes. The residues that fluctuate and diverge from the overall structure are represented by peaks in the RMSF plot (Fig. 4). Interestingly, it was seen that RMSF values of all FTO-ligand complexes exhibited an overall lower or similar RMSF value when compared to apo FTO protein during the simulation time of 100 ns. The average RMSF values observed were 0.18, 0.14, 0.18, 0.16, 0.16, and 0.19 nm for FTO-protein, FTO-BT012 complex, FTO-BD064 complex, FTO-RL442 complex, FTO-RL074 complex, and FTO-DTcomplex, respectively.

The compactness and rigidity of the FTO protein and FTO-ligand complexes were investigated using, Rg (Fig. 5). The average Rg value of FTO protein, FTO-BT012 complex, FTO-BD064 complex, FTO-RL442 complex, FTO-RL074 complex and FTO-DT complex were noted to be 2.54, 2.59, 2.54, 2.53, 2.52, and 2.57 nm, respectively. The result from the Rg profile showed that FTO ligand complexes along with the reference remained stable between 2.5 to 2.65 nm throughout the 100 ns time scale. It was observed that FTO-RL442, FTO-RL074 complex and FTO-BT012 complex showed lesser gyration and closed confirmation when compared to other complexes.

The formation H-bond is the critical driving force that determines the stability and drug specificity. It was observed that from MD trajectories results of H-bond analysis, fluctuation of H-bonds was between 0 and 4 for FTO-BT012 complex, 0 and 4 for FTO-BD064 complex, 0 and 3 for FTO-RL442 complex, 0 and 5 for FTO–RL074 complex and 0 and 2 for FTO-DT complex (Fig. 6). The average number of H-bonds in brown algal complexes FTO-BT012 and FTO–BD064 were 1.13 and 0.36, respectively, and these complexes had formed a maximum of four bonds. Similarly, the average number of H-bonds in red algal complexes FTO-RL442 and RL074 were 0.20 and 1.73, respectively, in addition FTO-RL074 complex had reached a maximum of five H-bonds.

Further to verify the post simulation interaction study, three best complexes FTO-BT012, FTO-RL074 and FTO-BT012 were selected along with the control FTO-DT complex (Fig. 7). The final time frame of 100 ns structure of these complexes were dumped and analysed based on the H-bond interaction pattern. The post simulation analysis of FTO-BT012 and FTO-RL442 complex reveals that H-bond between Ser229 is present at the binding site of FTO whereas FTO-RL074 reveals H-bond between Thr92 and Arg96. The residues Thr92, Arg96, and Ser229 are found to be present in the active site of FTO inhibitor (DT).

Binding free energy calculations

To understand the subsequent binding of ligand-protein complexes from dynamic simulation, binding free energy was calculated using the MM-PBSA approach. This method calculates the binding free energy by adding all energy types such as electrostatic energy, polar solvation energy, SASA (Solvent Accessible Surface Area) energy and van der Waals’ energy (Table 2).

The binding energies for FTO-BT012 complex, FTO-BD064 complex, FTO-RL442 complex, FTO-RL074 complex, and FTO-DT complexes were subsequently measured for last stable 20 ns MD trajectories. The total binding free energy of the reference molecule FTO-DT complex was −46.24 kJ/mol. Among other complexes, FTO-BT012 exhibited lowest binding energy of −102.06 kJ/mol whereas FTO-BD064 complex, FTO-RL442 complex, and FTO-RL074 complex possessed binding energy of −61.27, −57.40, and −80.24 kJ/mol, respectively. The binding free energy of FTO-ligand complexes was lower than DT. FTO-RL074 possessed lower van der Waals’ interaction of −105.94 kJ/mol in comparison with reference molecule that provides stable binding. Also, the FTO-RL074 complex exhibited highest SASA energy value of −13.35 kJ/mol in comparison with DT. In addition to this, the FTO-RL074 complex showed highest electrostatic energy of −86.45 kJ/mol that played a vital role in binding of protein-ligand complex.

The residues that contributed to the overall binding are alone depicted in the Fig. 8 for clear representation of the values. The residues that were majorly involved in energy contribution for all complexes were Val83, Ile85, Leu90, Thr92, Pro93, Leu109, Tyr220 and Met226 (Fig. 8). Interestingly, the abovementioned residues contributed to the active site of protein that plays a major role in maintaining the stability of the FTO-ligand complexes.