In silico molecular docking and molecular dynamic simulation of agarwood compounds with molecular targets of Alzheimer’s disease 

Selection and Preparation of ligands

Phytocompounds from the Aquilaria plant species were selected based on the previous literature and their structures were retrieved from the PubChem database (refer Figure 1). The phytocompounds structures were prepared by adding polar hydrogens, Gasteiger charges and by performing energy minimization in UCSF Chimera 1.16 using default parameters.³⁶

Figure 1. 2D structures of agarwood phytocompounds chosen.

Prediction of physicochemical properties

Prediction of pharmacokinetic and pharmacodynamic features can be accomplished by using the physicochemical characteristics of chemical substances, such as lipophilicity (LogP), solubility (LogS), and polar surface area and volume (PSA).³⁷ It is crucial to analyze these features because they affect how they interact with transport proteins and enzymes that are involved in drug clearance.³⁷ In the present study, we have used Molsoft tools (https://www.molsoft.com/mprop/) to predict the number of H-bond donors (No. of HBD) and acceptors (No. of HBA) present, polar surface area (MolPSA), lipophilicity (Mol LogP), solubility (Mol LogP), and the molecular polar surface area and volume (Mol PSA) of tested ligands.³⁸ Also, ADMET properties such as absorption, digestion, metabolism, excretion, and toxicity properties of selected ligands were predicted using ADMETlab 2.0 (https://admetmesh.scbdd.com/).

Selection and Preparation of receptors

The crystal structures of different proteins implicated in the pathogenesis of AD were retrieved from RCSB PDB structural database (https://www.rcsb.org/). The stereo-chemical properties, Ramachandran graph and values of selected proteins were evaluated by Molprobity server.³⁹ Chimera 1.16 (RRID:SCR_004097) was used to generate any missing residues in the selected target proteins. Following the removal of unneeded nonstandard heteroatoms, polar hydrogens and Gasteiger charges were added. All targets' structural details were refined using the steepest descent and conjugate gradient algorithms (100 steps each) with amber force field (Amber ff14SB).⁴⁰ Then, using AutoDock tools 1.5.7 (RRID:SCR_012746), the energy-minimized protein structures were transformed into ‘pdbqt’ format. A list of proteins along with their PDB IDs are given in Table 1.

Protein-ligand docking

Docking was performed with Autodock Vina as described in our previous study.⁴¹ The grid box's dimensions were fixed at $\mathit{XYZ}=30\AA \times 30\AA \times 30\AA$ which was found to be the best size for the default exhaustiveness (=8), and the ligand binding site was positioned in the middle of the grid box. The spatial dimension (XYZ axis) and the grid box’s size were specified in a configuration file. Using AutoDock vina version 1.1.2's (RRID:SCR_011958) command line interface, docking was accomplished. The obtained results are restricted to nine binding modes. The log file created included a list of the increasing binding modes and their associated binding energies. The BIOVIA Discovery studio visualizer 2021 was used to view the binding modes. All non-bonded interactions were recorded (DOI: dx.doi.org/10.17504/protocols.io.3byl4j362lo5/v1).

Molecular dynamic (MD) simulation

MD simulation using Gromacs 2020.5 (RRID:SCR_014565) was used to track the structural stability of the docked complexes. Gromos96 force field was used to create the topology of the protein, and PRODRG server (http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg/submit.html) was utilized to create the topology of the ligand. The docked complexes were solvated in a cuboidal box with adequate size to fit the complete complex in the middle using a “Simple Point Charge“ (SPC) water model. Appropriate number of counter ions (Na+/Cl-) were used to neutralize the simulated systems. The undesirable contacts and steric conflicts were then removed from the neutralized systems using steepest descent followed by conjugate gradient methods for 50,000 steps each.

The NVT ensemble used to maintain constant number of atoms, volume, and temperature, further NPT ensemble was used to maintain constant pressure. In this study, we set temperature and pressure constant at 300K and 1 bar respectively. Further, followed by 1ns of equilibration, unrestrained MD simulation was performed for a period of 100ns in solvent. The Particle Mesh Ewald (PME) method was used to handle coulomb electrostatic interactions, while the LINear Constraint Solver (LINCS) algorithm was used to limit H-bonds. Using a cut-off value of 14 Å, the non-bonded contacts were trimmed. The trajectories generated were analyzed using some of the inbuild gromacs tools like ‘gmx rms’, ‘gmx rmsf’, ‘gmx hbond’, ‘gmx gyrate’, ‘gmx sasa’, etc. and other additional packages for specific analysis wherever required. Conformational changes at the secondary structural level were monitored by using Dictionary of Protein Secondary Structure (DSSP) software (RRID:SCR_002725) (DOI: dx.doi.org/10.17504/protocols.io.36wgqjmp5vk5/v1)

Molecular mechanic/Poisson-Boltzmann surface area (MM-PBSA) calculation

MM-PBSA in conjunction with MD simulations is commonly used to determine the binding free energy of protein and ligand complexes.⁴² It uses the following equation as:

Enrichment analysis of protein targets

In STRING database (version 11.0), the proteins Glycogen Synthase Kinase 3 Beta (GSK3B), Monoamine oxidase-A (MAOA), Butyrylcholinesterase (BChE), Acetylcholinesterase (AChE), beta-site amyloid precursor protein cleaving enzyme 1 (BACE1), Legumain (LGMN), and Glutamate Ionotropic Receptor NMDA type subunit 1 (GRIN1) were searched for Homo sapiens. Only these proteins were selected for enrichment analysis, because most of the tested ligands showed good binding affinity only towards them. Gene ontology (GO) analysis was used to pinpoint the biological processes. Based on the available literature, the pathways contributing to AD pathogenesis were selected.

Network pharmacology

Cytoscape tool (version 3.9.1) was used to build the integrated network between chemicals, protein targets, and regulated processes. In the network analyser, the entire network was treated as “Direct” during network construction. The resulting network was inspected using a topological parameter “edge count.⁴⁴” The node size and colour were set using low values to small sizes and low values to bright colours, respectively.⁴⁵

Prediction of ADMET analysis

The physicochemical properties of all the chosen agarwood compounds were studied to gain more structural features of the individual phytocompounds (Table 2). Also, pharmacodynamic and pharmacokinetic properties of these compounds were analyzed using ADMETLab2.0. The analysis reveals that all the chosen compounds were less/not toxic. Additionally, it is interesting to note that most of the compounds had the capacity to cross the blood brain barrier (BBB). Though three phytocompounds namely, aquilarisin, aquilarisinin, and aquilarixanthone failed the Lipinski rule, they were shown to cross the blood brain barrier effectively and express lesser or no toxicity (Supplementary Material ADMET properties of agarwood compounds.xlsx). The blood brain barrier penetration values for agarotetrol, aquilarisin, aquilarisinin, aquilarixanthone, and pillion were predicted to be 0.8, 0.029, 0.35, 0.011, and 0.01, respectively. All the chosen ligands were predicted as non-carcinogenic and showed drug-likeliness.

Molecular docking analysis

In this study, docking of 41 agarwood compounds was done against 13 target proteins of AD. The selection of specific AD targets has been made considering their close and direct association with AD pathogenesis and its disease progression. Figure 2 represents heatmap of the binding energy estimated from the docking pose of agarwood compounds towards the tested molecular targets of AD. The scale of heat map ranges from least (blue) to highest binding (red) affinity was predicted based on the docking results. The binding interactions of only those ligands with least binding energy than the respective native ligands of AD target proteins are highlighted by circle in the heatmap (refer Figure 2). Ligands such as aquilarisin, aquilarisinin, and aquilarixanthone showed good binding affinity with most of the selected AD related proteins. Aquilarisin showied good binding affinity towards 12 of the 13 tested AD related proteins except with 4AU8. Aquilarisinin exhibited the highest binding affinity towards 8 of the 13 selected AD proteins including 1PBQ, 1Q5K, 1UDT, 4BDS, 4KKE, 4M0F, 5BTR, and 5LUA. On the other hand, aquilarixanthone was found to show highest binding affinity towards 9 of the 13 selected AD proteins including 1Q5K, 1UDT, 4AU8, 4BDS, 4KKE, 4M0F, 5BTR, 5LUA, and 6GZM.

Figure 2. The heatmap represents the binding energies (kcal/mol) of the agarwood compounds docked with target protein of AD.

Ligands with higher binding energies than the respective native ligands of target proteins are denoted in red circles.

From the heatmap it was clear that proteins such as 1Q5K (GSK3beta), 2Z5X (MAO-B), 4BDS (BChE), 4M0F (AChE), 5BTR, 5LUA (AEP) formed best docked complexes with most of the selected agarwood compounds. The site-specific non-bonded interactions of some of the representative top docked complexes showing highest binding affinity and conserved binding pocket interactions are listed in Table 3 (and also the 3D and 2D structures of top docked complexes are shown in Figure 3). Further, to gain more detailed insights to their structural stability and intermolecular interactions we selected some of the best docked representative complexes for MD simulations. The phytocompounds bound to these selected targets had been chosen considering binding energy, H-bonds and other nonbonded interactions including hydrophobic and electrostatic interactions. Additionally conserved binding pocket interactions had also been considered and compared with the respective known control inhibitors for the respective targets. Thus, we selected total seven complexes namely 1Q5K-AXN (complex1), 2Z5X-PLN (complex2), 4BDS-ANN (Complex 3), 4M0F-AXN (Complex 4), 5IE1-ASN (Complex 5), 5LUA-AGT (Complex 6), and 1PBQ-ASN (Complex 7) for MD simulations. It was observed that the agarwood phytocompounds bind to the AD targets to form stable complexes, and these complexes are maintained through H-bonds and other non-bonded interactions (as shown in Table 3).

Figure 3. Binding mode of the Agarwood phytocompounds to the AD targets (shown in 3D) predicted by docking study and their intermolecular interactions stabilizing the docked complexes in 2D.

Structural stability of the best docked complexes using MD

We selected 7 complexes of agarwood compounds 1Q5K-AXN (complex1), 2Z5X-PLN (complex2), 4BDS-ANN (Complex 3), 4M0F-AXN (Complex 4), 5IE1-ASN (Complex 5), 5LUA-AGT (Complex 6), and 1PBQ-ASN (Complex 7) having least binding energy towards selected AD targets. MD simulation quality check was performed over all the 7 trajectories by plotting temperature, pressure, and potential energy. During the 100ns simulation, the temperature and pressure were held constant at 300 K and 1 bar, respectively. We observed less fluctuations in the potential energy of the simulated systems, suggesting the well equilibration of all the complexes during simulation. In order to better understand the structural stability, the root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), and solvent accessible surface area (SASA) were also measured.

Root mean square deviation (RMSD) analysis

The backbone RMSD values were plotted over the trajectories revealing the stable dynamics expressed by all the simulated complexes during the 100ns. Figure 4 represents the backbone RMSD values of GSk3beta (indigo), BChE (Green), AChE (crimson), BACE1 (turquoise), and AEP (olive) showing moderate fluctuations up to 30ns (equilibration period), however, RMSD values are well stabilized after 50ns for these five proteins followed by equilibration period of 30ns. However, the RMSD of backbone atoms of NMDA (orange) and MAO-B (maroon) reached equilibrium after 60ns and 80ns during simulation. The average backbone RMSD values of GSk3beta, BChE, AChE, BACE1, AEP are observed as 2.5Å, 3.2Å, 2.1Å, 2.5Å, 2.4Å, etc. On the other hand, average RMSD values of NMDA, and MAO-B are found to be ~3.2Å and ~3.6Å, respectively. The representative MD simulation end structure extracted from the trajectory is compared with the initial starting MD simulation of the respective complexes. The analysis reveals the RMSD of structural superimposition of MD starting structure (0 ns) and end structure (100 ns) for studied 7 complexes has been found to be 1.339Å, 1.186Å, 1.122Å, 1.314Å 1.230Å, 1.161Å, 1.146Å. The ligand RMSD values also follows the similar trend (i.e., RMSD values <5Å) for all the ligands, however, the ligand7 (ASN) represents structural changes due to existence of torsions, which trigger the diverse ASN conformations during simulation. Thus, in general, all the 7 complexes express stable dynamics during 100ns MD simulation with backbone and ligand RMSD value of <5Å (refer Figure 4A and B).

Figure 4. RMSD of A) backbone and B) ligands of selected complexes.

Root mean square fluctuation (RMSF) Analysis

RMSF gives qualitative measure provides detailed insights to the conformational flexibility of protein structure. The RMSF value plotted for the C-alpha atoms of all the simulated systems is shown in Figure 5. The C-terminal region in all the complexes show highest fluctuations when compared to other regions in the protein structure. Careful observation of the three-dimensional structure reveals that increased RMS fluctuation values at C-terminal region are mainly due to the absence of native folded secondary structure. Overall binding pocket residues show the least RMS fluctuations (<3Å) as they actively participate in the stable non-bonded contacts. Also, other flexible loops and N-terminal regions express moderate to high fluctuations in the RMSF values. Due to a maximum residual fluctuation of up to 5, a greater peak between 175 and 190 amino acid residues was detected in complex 1. Only the C- and N-terminal portions of complex 2 exhibit persistent amino acid variation. Two higher peaks are observed in complex 3; one between 370 and 390 and one between 75 and 90 amino acid residues. Three higher peaks in complex 4 have been noticed at positions between 370 and 390 amino acid residues, between 250 and 275 amino acid residues, and between 40 and 60 amino acid residues (3.8Å, 3.9Å, and 3.8Å, respectively). Between 150 and 200 amino acid residues and 300 and 350 amino acid residues, complex 5 shows two higher peaks (each 5.5Å). A higher peak (4.5Å) has been seen in complex 6 between 90 and 110 amino acid residues. Three higher peaks (at 5Å, 6.5Å, and 4.5Å) have been discovered in complex 7, two of which were between 0 and 60 amino acid residues and one between 80 and 110.

Figure 5. RMSF of selected complexes.

Analysis of Radius of Gyration (Rg) and solvent accessible surface area (SASA)

The radiation of gyration (Rg) explains the protein folding/compactness of the molecule, hence we analyzed the compactness of the protein-ligand complexes and exposure of hydrophobic core of the protein to the solvent upon ligand binding. The variation in Rg and SASA of selected complexes is given in Figure 6. The Rg value for complexes 1-7 (except complex 2) is well stabilized while complex 2 shows steady decrease in the Rg value during the MD simulation. Also, SASA values of all complexes show significant structural stability and represent the formation of compact globular shape during 100 ns simulation. The average Rg and SASA values of complexes 1 to 7 range between ~1.75 to 2.5Å and 155 to 320 (nm2). Thus, we observed stable complex formation in all the complexes (refer Figure 6).

Figure 6. Compactness of docked complexes and exposure of hydrophobic pocket to the solvent estimated by Radius of Gyration (Rg) and solvent accessible surface area (SASA) respectively.

Intermolecular interactions observed in docked complexes

Formation of H-bonds between selected protein-ligand complexes has been monitored during 100 ns simulation. Figure 7 denotes the number of H-bonding interactions formed during simulation. Table 4 lists comparative analysis of non-bonded interactions observed in the initial starting structure (0ns) and MD simulation end structure (100ns). In general, the native non-bonded contacts were well conserved even after 100 ns MD simulation, revealing relatively more stable complex formation in all the simulated complexes. Hydrophobic interactions equally contribute in stabilizing these complexes. Moreover, the consistency of the observed H-bonds also supports the stable complex formation during simulation (Refer Figure 7). The snapshot of initial and final MD structure is shown in Figure 8, revealing the binding of ligands to the conserved binding site of their respective targets. Also, most of the complexes show compact folding during simulation forming much compact globular structure.

Figure 7. Number of H-bonds formed in 7 docked complexes during MD simulation.

Figure 8. The representative snapshot of the initial and final structure of the MD simulation.

Estimation of Binding Free Energy using MMPBSA

The binding free energy for all the simulated complexes has been quantitatively measured by using MMPBSA methods over the well equilibrated trajectory observed between 0 to 100 ns. Table 5 represents the energy components including molecular mechanics, van der Waal (vdW) interactions, electrostatic, polar and nonpolar energies that significantly contribute to the binding free energy. The estimated binding free energy for complexes 1 to 7 is -81.018 ± 61.364, -159.438 ± 10.190, -227.959 ± 13.745, -152.764 ± 15.897, -149.090 ± 16.646, -79.236 ± 19.623, and -146.796 ± 12.694 kJ/mol, respectively. It has been observed that binding free energy of protein-ligand complexes is significantly influenced by both electrostatic and vdW interactions.

Conformational changes at secondary structural level

We further examined the secondary structural changes during simulation period using DSSP. The complexes 1-7 using DSSP plot have been shown as panel A-G in Figure 9. Various components of secondary structures are shown in specific colors as shown in figure legend. We noticed that major secondary structural components such as alpha helices and beta sheets are relatively much stable and express less variations in the secondary structure. Also, interestingly, the architecture of binding pocket is well maintained throughout the simulation due to rigidity provided by the alpha helices and beta sheets in all the simulated complexes. Thus, overall binding of ligand does not affect the secondary structure of protein. However, minor variations such as shortening, or extensions were noticed at flexible loop and linkers connecting the alpha helices and sheets. The H-bonds provide structural stability to the alpha helices and beta sheets in all the complexes.

Figure 9. The dictionary of protein secondary structure elements.

A) complex 1, b) complex 2, c) complex 3, d) complex 4, e) complex 5, f) complex 6, g) complex 7.

Enrichment analysis and Network construction

The GSK3B, MAOA, BCHE, ACHE, BACE1, LGMN, and GRIN1 all influenced 18 different biological processes. Among them, 15 were associated with AD. Among 7 targets, 6 targets viz., BCHE, ACHE, BACE1, GSK3B, GRIN1, and LGMN were enriched for modulation of chemical synaptic transmission (GO:0050804) and 3 targets viz., BACE1, GSK3B, LGMN were enriched for cellular response to amyloid-beta (GO:1904646). Following these, regulation of neurotransmitter levels, neuron death, cell communication, signaling, synaptic plasticity, neuron apoptotic process, learning, etc were also associated with SK3B, MAOA, BCHE, ACHE, BACE1, LGMN, and GRIN1 targets. Table 6 represents the biological processes modulated by best docked phytocompounds. Among the 7 targets, GRIN1, GSK3B, BACE1, and LGMN scored the highest edge count within the network of and were involved in multiple biological processes for the regulation of AD (Figure 10).

Figure 10. Network representation of compounds, targets, and modulated molecular processes involved in AD.