In silico molecular docking and molecular dynamics examination of Andrographis paniculata compounds of Andrographolide, Neoandrographolide, and 5-hydroxy-7,8,2’,3’-tetramethoxyflavone inhibition activity to SARS-CoV-2 main protease

In this work, Andrographis paniculata compounds of Andrographolide, Neoandrographolide, and 5-hydroxy-tetramethoxyflavone inhibition activity to SARS CoV-2 main protease were examined through in silico molecular docking and molecular dynamics simulation, with Remdesivir as control ligand. Docking score and MMGBSA were examined as well as molecular dynamics parameters: RMSD, RMSF and Protein ligand contact fraction. Our study found that Andrographis paniculata compounds of Andrographolide, Neoandrographolide, and 5-hydroxy-tetramethoxyflavone have comparable inhibition activity to SARS CoV-2 main protease in comparison to Remdesivir. 5-hydroxy-tetramethoxyflavone has the lowest docking score, which was further validated by protein ligand contact fraction examination, although MMGBSA score is lowest for Remdesivir.

Computer aided drug designing approach, molecular docking, and molecular dynamics, has been very effective in the past in drug discovery for various diseases [14]. Through in silico molecular docking, Rajagopal et al., found that Andrographolide has significant bond with the active site of SARS CoV-2 main protease, PDB ID 5R82 and might produce significant activity [15]. Hiremath, et al. [16] also showed some compounds in Andrographis paniculata have potential as SARS-CoV-2 inhibitor. Main protease (Mpro) is viral protease that perform cleaving two polyprotein into individual non-structure protein. Mpro has important role in viral replication, thus can be target protein for inhibition [17]. Remdesivir has been reported to be a potential inhibitor of SARS-COV-2 Mpro and reported to be effective in treating COVID-19 [18][19][20]. In this paper, we investigate the inhibition activity of some Andrographis paniculata compounds (Andrographolide, Neoandrographolide, and 5-hydroxy--tetrametho-xyflavone) to Mpro of SARS-CoV-2 through molecular docking, as well as examine dynamics interaction between compounds and Mpro through molecular dynamics simulation. Molecular docking and dynamic simulation of Remdesivir with Mpro was also performed to be positive control.

Ligand and protein preparation
3D structure of sambiloto compounds and remdesivir as control, called as ligands, were obtained from PubChem database in the form of a 3D conformer with the .SDF format for later preparation using the Ligprep tool available on the Maestro Schrodinger software. The use of remdesivir as control ligand was determined as the study by McCreary et al. showed that some of the randomized clinical trials suggested that remdesivir may improve recovery for individuals hospitalized with COVID-19, although inconclusive findings exist [21][22][23]. As for the preparation process, the ligands were desalted, and hydrogen atoms were added. All possible ionization states and tautomer were created using the Epic configuration at a target pH of 7.0 +/-2.0. The lowest energy conformation is obtained by minimizing energy using OPLS 2005.
The structure of the SARS-Cov-2 (Mpro) protein was obtained from the protein data bank (PDB ID: 7AOL) in .pdb format [24]. The typical structure of the proteins in this .pdb file cannot be used directly in in silico docking simulation calculations because they might consist heavy atoms and may include native ligands, water molecules, metal ions and crystallized cofactors. The proteins were prepared using the Protein Preparation Wizard in the Maestro Schrodinger software by adding missing hydrogen atoms, side chains, and loops and removing water molecules from the target. Next step is adjusting the ionization and tautomeric state of the hetero group using Epik. The final step is optimization of the hydrogen bonding and refining the structure by limiting the RMSD minimization to 0.3Å.

Grid generation
Grid box in simulation was prepared using the Receptor Grid Generation tool. This was done to specify that the simulation process was carried out in a certain area. The grid is set at the centre of the inhibitor RQH-403 which is the native ligand of the PDB file 7AOL at the protein active site. The Van der Waals radius scale factor was set at 1.0 with a partial charge cut off 0.25. Then, the grid area is the centroid of workspace ligand, and the box size is 20 Å.

Molecular docking
Molecular docking was done using the Ligand Docking tool by loading the receptor file (.zip format) and the prepared ligand (.maegz format) into working space. The Van der Waals radius scale factor was set at 1.0 with a partial charge cut off 0.25. Extra precision (XP) algorithm and flexible ligand sampling was used. The output of this docking process will be saved in .mae format.

Molecular Mechanics-Generalize Born Surface Area (MM-GBSA)
The results obtained from molecular docking did not show the affinity or free energy of the bond between the protein and the ligand. Therefore, it is necessary to calculate the MM-GBSA to determine the binding free energy. MM-GBSA can be calculated using Prime MM-GBSA tools on Maestro Schrodinger by loading the output file of molecular anchoring (.mae format). The solvation model in the option was set as VSGB and the sampling method on protein flexibility was minimized.

Molecular Dynamics
Prior to molecular dynamics simulation, it is necessary to prepare a system to determine parameter settings using the System Builder tool. Solvent model in the system was set as simple point charge (SPC) model with ensemble NPT. The simulation box shape was set to orthorhombic with distance a = 10.0 Å b = 10.0 Å c = 10.0 Å with minimized volume. Ion neutralization was carried out by adding Na+ ions. Salt was added with concentration of 0.15 M.
Molecular dynamics simulation was set for 10 ns to all ligands with the trajectory interval being of. Ensemble class uses NPT with a constant temperature of 310 K and a constant pressure of 1.0132 bar. The output obtained from this simulation was visualized for discussion.

Molecular docking and MM-GBSA
In silico molecular docking involves conformation and orientation of the ligand when it binds the target. The docking score is a quantitative description of the bonding modes that occur in molecular binding and the score can distinguish whether ligands bind the target or not [25]. The lower the energy value indicates that the intermolecular bonds forces are stronger.
To add discussion, we incorporated MM-GBSA (molecular mechanics-generalized Born and surface area) or binding energy, which is combination of the probability of free ligand-protein and the effect of the solution, composed of coulomb energy, covalent binding energy, Van der Waals energy, lipophilic energy, generalized Born electrostatic solvation energy, hydrogen binding -According to Table 1, most of the compounds contained in the sambiloto have lower docking score than remdesivir (-1.204 kcal/mol). 5-Hydroxy--tetramethoxyflavone compound was observed to has lowest docking score at -3.180 kcal/mol, while the lowest binding energy is Remdesivir at -45.99 kcal/mol. Difference between docking score and MMGBSA due to the different approach where MMGBSA allows for complex free energy decomposition into contributions from various groups of atoms and types of interaction [26]. The ligands-target interactions were shown in Figure  1-4 for Remdesivir, Andrographolide, Neoandrographolide, and 5-hydroxy--tetramethoxyflavone respectively.
For Remdesivir-target interaction, Glu 178 is an acceptor that forms two hydrogen bonds with the hydroxyl and amine groups of Remdesivir (Fig.1).For Andrographolide-target interaction, they form two hydrogen bonds. One OH group is a donor to Phe 103, and another is acceptor with Lys 102 (Fig. 2).For Neoandrographolide-target interaction, they form three hydrogen bonds to target. donor to Gln 110, donor to Arg 105, and acceptor to Phe 103 (Fig. 3).

Fig. 2. 2D projection of Andrographolide interaction with Mpro
While, for 5-Hydroxy--tetramethoxyflavone ligand interaction, Arg 105 interacts with two groups at once (Fig. 4). It Formed hydrogen bonding and Pi-cation and the Glu 178 also binds to two hydroxy groups in ligand and acts as an acceptor.  Neoandrographolide has the lowest score since the net of exchange energies between two residues (Arg 105 and Phe 103) in favour of stronger interaction.

Fig. 4. 2D projection of 5-Hydroxy--tetramethoxy flavone interaction with Mpro
The difference between docking score and MMGBSA values exist, and it is not supposed to be al. [27], found there are almost no correlation between those two scoring functions, and both can be used independently as scoring function to predict binding affinity and ligands performance. According to our results, all of three sambiloto ligands, has comparable performance with Remdesivir, and in the case of docking score, all sambiloto ligands outperform Remdesivir.

Molecular Dynamics (MD)
Now we discuss about molecular dynamics simulation interaction between target and ligands. Molecular dynamics simulations were carried out to determine the stability of the ligands when interacting with proteins by considering the effects of the applied ensemble [28]. The results of the molecular dynamics simulation are the root mean square deviation (RMSD), root mean square fluctuation (RMSF) and protein-ligand contact [29].
RMSD shows whether the system has reached stability by observing the fluctuation until the value is not significant enough after reaching a specific period. If the RMSD is still fluctuating significantly until a specific time, then the system is not stable. The RMSD of the C and the side chain of the protein along simulation with sambiloto ligands were shown in Fig. 5-8. RMSD value is used is used to compare the docked conformation to reference conformation. The RMSD values below 2.0 Å with respect to the corresponding atomic structures are considered favourable [30]. In the beginning of the simulation, there was a sharp increase in the RMSD value for all compounds where the interaction attempted to reach equilibration. To reach a stable condition (equilibration), remdesivir took about 1.5 ns and during simulation time it reach up 2 Å, but only a glimpse. while andrographolide is less than 1 ns. In neo-andrographolide, stability occurred after 2 ns. while for the 5-hydroxy--tetramethoxyflavone compound in sambiloto, the RMSD of C Å and eventually reach stability under 2 Å after took 8.5 ns to reach stability. But none all the compounds have RMSD value of side chain below 2.0 Å while RMSD value of C were below 2.0 Å for all ligands. We had extended the time, but we could not find interesting features    The stability of the ligand-protein interaction can be examined too through the RMSF values. In contrast to RMSD, RMSF is calculated for each residue of to see how much fluctuation of each residue during the simulation [31]. In general, RMSF provides fluctuations of individual atoms throughout the simulation [29]. RMSF of all ligands-protein simulation were shown in Fig. 9-12, with green line indicates interacting residue.    Based on the RMSF evaluation, some of the protein residues of remdesivir and sambiloto MD simulation showed fluctuations higher than 2.5 A. Further observation shows that the residue fluctuated more than 2,5 Å were not the interacting residue. Whereas all interacting residues, indicated by green line, generated RMSF below 2.5 Å, thus it is considered acceptable [32] indicating that the interaction of ligands with proteins are stable.

Protein-ligand contact fraction
In this subsection, protein-ligand contact was investigated to further validate docking scores and MM-GBSA values. Protein-Ligand contact during the simulation was recorded and processed as a simulation interaction diagram and timeline representation. These results inform the type, fraction, which residue is involved. Interaction fraction value inform fraction of simulation time that a specific interaction occurs. For instance, interaction fraction 0,6 means that a residue interacted with ligand for 60% of simulation time, or in this study of 10ns molecular dynamic equal to 6ns. Value of more than 1,0 is possible for a single residue as it may contact multiple groups in a ligand.    Most of ligands interact with Arg 105 residues through various interactions except for Neoandrographolide with its dominant interaction with Arg 298 residue. However, it is for 5-hydroxy--tetramethoxyflavone which has high frequency of hydrophobic interactions, distributed to Arg 105, Phe 134, and pro 184 residues, these hydrophobic interactions which distributed, are more frequent in high-efficiency ligands [33], thus validates the lowest docking score of 5-hydroxy-tetramethoxy flavone.

Conclusion
In this study, we have conducted in silico molecular docking and molecular dynamics simulation of Remdesivir (as control ligand) and three sambiloto compounds: Andrographolide, Neoandrographolide, 5hydroxy--tetra-methoxyflavone. We found that all the investigated sambiloto compounds have the lower docking score than Remdesivir, while 5-hydroxy--tetramethoxyflavone is the lowest. although MMGBSA results showed that Remdesivir still has lowest score among all examined ligands.
We also found that RMSD value of C for all ligands consistently below 2.0 Å for all ligands, although for RMSD value of side chain for all ligands are still higher than 2.0 Å. However, RMSF results demonstrated that all the ligand's interactions are dynamically favourable where the RMSF of interacting residues are below 2.5 Å. The protein ligand contacts fraction examination validated the lowest docking score of 5-hydroxy-tetramethoxyflavone.