Hepatitis B is a viral infection that attacks the liver and can cause both acute and chronic disease. WHO estimated that 257 million people are living with hepatitis B virus infected, which is defined as hepatitis B surface antigen (HBsAg) positive. Hepatitis B caused about 887.000 deaths, mostly from complications, including cirrhosis and hepatocellular carcinoma in 2015.1 Chronic hepatitis B infection remained a significant threat to public health despite the fact HBV vaccination programs have been extensively implemented in the past. The approved agents of chronic HBV infection treatment are interferon and nucleoside analogs. However, the current treatment still less effective because of the side effect of interferon and the viral resistance of nucleoside analogs. So, it is necessary to find and develop new drugs with novel structures and mechanism.2-3 It is therefore important to continue research to identify new anti-HBV targets and discover novel compounds with a unique mode of action for controlling viral replication as well as preventing drug resistance and its complications in the long term. One approach is to target virus assembly that was formed by dimerization of core protein.3-10
HBV is a small virus whose genome has four open reading frames which correlate with a complexity of function for viral protein. HBV core protein contains 183 residues that self-assembles to form the viral capsid. In an infected cell, it modulates almost every step of the viral replication process. Core protein bound to nuclear viral DNA and affect its epigenetics correlates with RNA specificity, assembles on reverse transcriptase-viral RNA complex specifically, participates in the regulation of reverse transcription, signals completion of reverse transcription to support virus secretion, carries both nuclear localization signals and HBsAg binding site.10 The HBV core protein is no known related to proteins in human cells. It is an excellent target for the development of new, virus-selective, safe and effective antiviral agents to improve treatment option for hepatitis B disease. It is consists of 183-185 amino acid that form N-terminal for capsid assembly (amino acid 1-149) and C-terminal domain for nucleic acid binding (amino acid 150-185). The viral capsid infectious HBV particles are formed from 120 copies of assembled core protein dimers enclosing the viral DNA. Small molecules which can disrupt functional HBV capsid assembly can be potent replication inhibitors of HBV.11
In most cases, the initial lead compound was found by process of random screening. Many classes of non-nucleoside compounds have been shown promising activity against HBV. In the last several years, in particular, the compounds obtained from natural sources prepared by synthesis/semi-synthesis have been reported as novel HBV inhibitors.2 One method used in the screening process is to use structure-based search is molecular docking. The main purpose of molecular docking is to understand the mechanism of drug works at a molecular level by observing drug interaction with the receptor. Computational approaches in docking simulation of drugs or candidate drugs into macromolecules as a target receptor then scoring of their potential complementary to the binding site are widely used in hit identification and lead optimization.12
Emodin (1,3,8-trihydroxy-6-methyl-9,10-anthraquinone) is derived from herbal medicines. It has been derived such as from extract of Phragmanthera austroarabica,13 Polygonum multiflorum,14 Cassia occidentalis,15 and Bryonopsis laciniosa fruit.16 Emodin has been reported as a weak activity but a persistent inhibitory effect on HBV replication, both in vitro and in vivo evaluation. It inhibits the secretion of HBsAg and HBeAg in the cultured cells. It inhibited the production of HBV DNA dose-dependently with IC50 was 0.21 g/L and the selected index CC50/IC50 was 1.95.17 Emodin and Astragalus polysaccharide (APS) reduced serum HBV DNA content of transgenic mouse models significantly, although this effect was weaker than lamivudine group, but lasted longer.18
Due to its activity, emodin is chosen as a lead compound for developing a new hepatitis B drug candidate because its structure showed a useful pharmacological activity and can be as the starting point for drug design. Furthermore, it had a high affinity for phospholipid membranes and could pass the cellular membranes easily because it has high phospholipid/water partition coefficient.19
By modification structure where one or more particular functional group of the molecule is removed or altered, it is possible to find out the essential groups. This study aims to explore the potential of emodin derivatives in inhibition of HBV core protein by molecular docking and dynamic simulation. Modification of emodin was designed by substitution a benzoyl, ortho-, meta, and para-toluyl group at the C3 position which replaces the hydroxyl group. Because the 1,8-hydroxy groups of Emodin are adjacent to the 9-carbonyl group, two intramolecular hydrogen bonds are thereby formed, while the 3-hydroxy group is far from the 9-and 10-carbonyl group, the 3-hydroxy is easier to convert than 1- and 8-hydroxy moieties. Of its results expected to be known later to be developed as potential candidate compounds inhibiting HBV replication with the target protein core. Therefore, this study conducted a study to look at the molecular docking and molecular dynamics of emodin compounds that have been modified by adding a benzoyl (FE-3), o-toluyl (FE-4), or p-toluyl(FE-6) group at C3 position replaces the hydroxy group.
Protein structure data and preparation
The three dimensional (3D) crystal structure of core protein hepatitis B virus N-terminal domain bound to NVR-010-001-E2 (PDB ID: 5E0I) was retrieved from the Research Collaboratory for Structural Bioinformatics (RCSB) Protein Data Bank (http://www.rcsb.org/). The crystal structure has a resolution of 1.95 A. It also has a structural weight of 109089.31 Da and amino acid length of 157 and contains six chains (Chain A, B, C, D, E, and F). The protein was imported to Molegro Virtual Docker (MVD) program. For molecular docking purpose, all the water molecules were removed because they were not considered during the scoring.
The 2D structures of emodin (1), 3-benzoyl emodin (FE-3), 3-ortho toluyl emodin (FE-4), and 3-para-toluyl emodin (FE-6) were drawn in Marvin Sketch. The structures were converted to 3D format as the most stable conformation and saved as Sybyl.mol2 file format for docking purposes.
The molecular docking program was performed using Molegro Virtual Docker (MVD) which incorporates highly efficient PLP (Piecewise Linear Potential), and MolDock scoring function provided a flexible docking platform. Validation of methods done by extracting ligands already present in proteins (code: 5J6) and performed the docking process. Programs that can return poses below RMSD value of 1.5 or 2 A is considered to have performed successfully.20 Further ligands present are replaced with an emodin and benzoylated emodin by performing alignment and docking simulations performed for each of these compounds. There are six chains of HBV core protein in 5E0I which each charged with its ligand. Docking parameters were set to 0.30Å as grid resolution, maximum iteration of 1500 and maximum population size of 50. Energy minimization and hydrogen bonds were optimized after the docking. Simplex evolution was set at maximum steps of 300 with neighborhood distance factor of 1. Binding affinity and interactions of ligands with protein were evaluated from the internal ES (Internal electrostatic Interaction), internal hydrogen bond interactions and sp2-sp2 torsions. Post dock energy of the ligand-receptor complex was minimized using the Nelder-Mead Simplex Minimization (using non-grid force field and H-bond directionality). From the reranking score, the best interacting compound was selected from each dataset.
The crystal structure of complex and ligands poses were taken and extracted from the docking workspace. Both the force parameterization and hydrogen atom compensation were implemented using the LEaP module. Subsequently, the force parameters of NVR-010-001-E2, emodin and emodin benzoylated were chosen through the GAFF force field. Since the AM1-BCC (AM1 with bond charge corrections) was found to perform well, the atomic partial charges were assigned through the AM1-BCC fitting technique located in the sqm module using the Antechamber module. Each of the systems was soaked in a truncated octahedral periodic box of TIP3P water which extended a distance of 10 Å to the closest atom of the complex, and an appropriate number of ions were added to neutralize the system. First, every system was energetically minimized by a total of 200 steps with the first 50 being steepest descent. For implicit solvent model, it was used the GB model of Hawkins, Cramer and Truhlar.22-23 Next, the solvated complexes were equilibrated by heating and density equilibration with weak restrain during 50 ps respectively then followed by 500 ps of constant pressure equilibration at 300K. All simulation would be run with shake on hydrogen atoms, a 2 fs time step and Langevin dynamic for temperature control. Finally, 20 ns of MD simulation was performed for the system. All the above operations were performed with the PMEMD module.22 System preparation, simulations, and following analysis were all carried out by using Amber12 package.
Amber analysis toolkit utilities were used to analyze MD trajectories produced during the last 20ns of the production run to analysis root mean square deviation (RMSD) and root mean square fluctuation(RMSF). During the simulations, 3D coordinate snapshots of protein-ligand complexes were extracted at 0.5, 5, 10, 15 and 20 ns that used to analyze the various interactions for each protein-ligand system using MVD. For RMSD and RMSF graphs plotting, were used MS Excel (2013). The number of distinct hydrogen bonds formed between specific amino acids residues and ligand atoms were determined to utilize the energy set at a maximum of -2.5 kcal/mol.
Docking studies and validation of scoring function
The RMSD values between the generated poses and the bound ligand of crystal structure 5E0I are given in Table 1. It was found that poses generated using Moldock score function of chain D show lower RMSD value of 0.59Å. The binding site was set with coordinates of x = 146.31; y = 6.43; z = 123.97 in radius/dimension of 15 Ǻ. The reference ligand NVR-010-001-E2 compound lay on the floor of the binding area and interacted with protein via hydrogen bond with the Trp 102 and Ser 106 residues.
Molecular dynamics simulation result analysis
The ligands interacted with amino acid residues with 2,5 kcal/mol were showed in Table 3.
The validity of the docking settings was assessed by the ability to reproduce the X-ray crystal bound conformation (pose). The co-crystallized ligands of NVR10-001E23D coordinates were used for all docking experiments. In ordered to validate the scoring function, before docking the screened molecules for selecting prospective hits, we have docked into 5E0I protein structure in MVD. There are 6 chains of HBV core protein in 5E0I were each charged with its ligand. Subsequently, we compared the conformation and position with the bound ligand conformation measured regarding the root-mean-square-deviation (RMSD). Emodin and its derivatives structures are docked to the receptor at the same place and coordinates previously. The docking scores and interactions between ligand-receptor were analyzed.
Interactions between emodin and derivatives indicated by hydrogen bonding and steric interaction between the ligand with the amino acid residues in the HBV core protein. Emodin interacted via hydrogen bonding with amino acid residue Phe 23 and Tyr 118 at chain D, and Ala 132 at chain E. FE-3 interacted with residue Ser 106, Trp 102, Pro 138, Tyr 118 at chain D and Ala 132 at chain E. FE-4 interacted with residue Leu 140, Tyr 118, Trp 102 and Ser 106. FE-6 interacted with residue Tyr 118.
The target interaction with emodin derivative ligands obtained from molecular docking is verified again by simulating molecular dynamics to see the stability of the interactions at a certain time. Dynamics simulation is a form of computer simulation where atoms and molecules interact over a period of time-based on the laws of physics. Newton’s equations require a short time step (femtosecond level) to be integrated to make the simulation short, generally using a nanosecond unit. Simulation of molecular dynamics is done using Amber program. All the systems in the study were well equilibrated and remained stable throughout the simulation time. Based on the simulation results of molecular dynamics, obtained RMSD and RMSF profiles, ligand-protein interactions and interaction changes that occur due to the treatment specified in the simulation.
It is seen that the RMSD profile initially increased drastically to a value of about 3 Å in a short time span, which is about 0-5 ns. After that, RMSD increases volatile, but the fluctuation value is not too large. At the end of the simulation of 20 ns, the RMSD value reached 5 - 6 Å. The system appears to be convergent, but overall there is still a fluctuated in RMSD value. The RMSD profiles of emodin derivatives and HBV core protein complexes could be shown at Figure 2.
The RMSF value (Root Mean Square Fluctuation) was evaluated to see fluctuations in each of the protein amino acid residues during the simulation, which illustrated the residual flexibility. Residues with high fluctuations show high flexibility, meaning interactions or unstable bonds. The root-mean-square fluctuation (RMSF) values of the systems were also calculated and plotted to compare the flexibility of each receptor amino acid residue. The RMSF profiles showed flexible regions of the systems during the simulation time. For all the systems, high fluctuations occur in the loop regions as compared to that in the other regions. A low RMSF value indicates the well-structured regions, while the high values indicate the loosely organized loop or terminal domains. Moreover, to investigate the motions of the important residues, interacting with the inhibitor in the binding site, the root mean square fluctuations (RMSF) for all important residues were calculated and depicted in Figure 3. This focused RMSF plot shows that none of the important residues present in the active site has an RMSF value more than 0.15 nm. The result has confirmed that the residues present in the active site were minimally distorted upon binding of the substrate, hit molecules and any of the known inhibitors. On the charts, visible flexibility varies on each residue. It is seen that the residues 73 and 140 have the highest fluctuations in all bonds with the emodin derivative ligands.
Hydrogen bonds are key players and directly responsible for affinity and specificity in the protein-ligand complexes. Thus, to evaluate the binding affinity of the ligands with the target protein, we calculated the hydrogen bonds between ligand and protein during 1, 5, 10, 15 and 20 ns of MD simulation. The graphic of the rerank score of emodin derivatives into HBV core protein during dynamic simulation was shown at Figure 4.
Based on the results of the dynamic simulation showed that the interaction between ligand with receptor changes and fluctuation. This is possible because the ligand position is in the intermediate regions of the dimers, rather than in the protein cavity so that the ligand is free to move even more with the presence of solvents and temperature rise. This still allows the occurrence of inhibition of viral replication due to a change in orientation of capsid assembly by core proteins.
As determined by docking studies, additional chemical groups are necessary for emodin to increase its potency, especially if they are designed to interact at the core protein dimer-dimer interface. FE-4 was predicted has potency as replication inhibitor according to docking and dynamic simulation Thus, some chemical modifications to address the identified undesirable properties may be necessary. We believe that our findings provide a suitable starting point for finding a new candidate anti-hepatitis B agent an following with in vitro and in vivo analyses.