Introduction

Hepatitis B virus (HBV) is a disease of global concern. This virus has infected approximately 400 million people worldwide, including 93 million people in China1,2. HBV is a highly variable virus with mutations that occur because of a shortage of proofreading activity of the DNA polymerase, which results in the formation of different genotypes and subtypes3. There are eight hepatitis B virus genotypes (A-H) throughout world. Genotypes B and C are the most common genotypes in Asia4. Compared to genotype B, genotype C is associated with later HBeAg seroconversion, more severe liver disease and a higher incidence of hepatocellular carcinoma (HCC)5,6,7.

Currently, the available antiviral drugs include interferon-alpha (IFN-α) and nucleotide analogue inhibitors of HBV polymerase, among which lamivudine (3TC), telbivudine (LdT), entecavir (ETV), adefovir (ADV) and tenofovir (TDV) have been approved by the FDA8,9 (Figure 1). ADV is an acyclic nucleotide analog that inhibits the HBV DNA polymerase and suppresses HBV replication at the level of DNA synthesis. It is efficacious even for the treatment of lamivudine-resistant chronic hepatitis B10,11.

Figure 1
figure 1

Chemical structures of five nucleotide analogue HBV polymerase inhibitor approved by US-FDA for treatment of HBV infection.

Full size image

To treat patients with chronic hepatitis B, a personalized approach allows for the treatment efficacy to be optimized, which involves individualizing the therapy according to the genetic differences of the infectious agents. Recently, increasing evidence has shown that there are differences in the antiviral response among HBV genotypes12,13,14,15,16. IFN-α therapy was reported to be less effective in European patients infected with genotype D compared to those infected with genotype A, whereas in the Asian population, genotype C showed worse performance than genotype B. A clinical trial of 3TC has demonstrated that patients infected with HBV genotype A had higher DNA levels than those infected with genotype D17. For ADV therapy, the data associated with the sensitivity among HBV genotypes are subtle, particularly between the two predominant genotypes in Asia, ie, B and C. Zeng et al identified a genotype-specific outcome, indicating that the HBV genotype B shows a greater reduction of viral DNA than genotype C (3.6 vs 3.1 logarithm copies/mL) at 48 weeks ADV post-treatment in a Chinese Han population with chronic HBV18. However, another study of ADV therapy reported that a difference in the level of viral DNA reduction (3.42 vs 3.65 logarithm copies/mL) was observed after 48 weeks of treatment, although the association with HBV genotypes was not statistically significant19. The differences in the average reduction of viral DNA could partially reflect the difficulty in analyzing the antiviral response to ADV therapy because the different ethnicities, endpoint definitions, sample size and whether the patients have been treated with other therapies would all affect the results of the clinical trials20.

Recently, a computational study of the possible mechanisms of the different inhibitory effects of benzothiadiazine molecules on HCV genotype-1a and genotype-1b NS5b polymerase was reported, in which the ligand sensitivity was investigated with homology modeling and molecular docking21. Clearly, the structural information of the HBV polymerase from genotypes B and C would also be helpful for investigating the molecular basis of ADV sensitivity. Because the crystal structure of the HBV polymerase is not available, homology models based on the X-ray structure of HIV-1 RT have been reported and utilized for determining the molecular mechanism of inhibition and drug resistance8,22,23,24,25. Thus, no computational study was reported regarding the molecular basis of drug efficacy among different HBV genotypes. In this study, the differences between genotypes B and C that influence the binding affinity of ADV are studied using computational approaches. Two 3D models were first constructed for the HBV polymerase from genotypes B and C. Molecular docking, dynamics simulations, and binding free-energy calculations were then performed to explore the binding patterns of ADV and provide insights into the difference in the binding free-energy of genotypes B and C.

Materials and methods

Representative sequence selection for different genotypes of HBV polymerase

HBV genotype classification is based on a nucleotide sequence divergence of 8%4, although in sequences of the same genotype, naturally occurring polymorphic sites are often observed. Therefore, consensus genotype reference sequences of HBV polymerase for genotypes B and C are essential for homology model construction. Using the NCBI virus sequence genotyping tool26,27, standard reference sequences of genotypes B and C were obtained. The resulting sequences were then analyzed using the online alignment tool, ClustalW228, to determine the representative sequences for the two genotypes.

Homology modeling of HBV polymerase-DNA complex

The genotype consensus reference sequences reported by Rhee et al9 were used to represent the HBV polymerase genotypes B and C. The sequence alignment between the two genotypes was performed using the EMBOSS online program29. Many HIV-1 RT crystal structures are available in the PDB database. Because the complexed ligand, TDV, is a methyl derivative of ADV, the complexed structure of HIV-RT/TDV (PDB ID: 1T05) was selected as the template30. The recently published alignment of HBV polymerase and HIV-1 RT that was developed by Daga et al25 was used as the reference template for modeling the HBV polymerase domain. Following the alignment, the homology models of genotypes B and C were built using the automated modeling module, MODELLER31, in the Discovery Studio (DS) software package version 3.0 (Accelrys, San Diego, CA, USA). The optimal level of the entire protein and loops was set to high and medium, respectively. The best structure for each genotype was then selected according to the PDF (probability density function) scores and the three-dimensional fit to the template protein. After the initial model was prepared, the coordinates of the template-primer DNA duplex (11/7), two Mg2+ cations and one original ligand [tenofovir diphosphate (TDV-DP)] were transferred from the template structure, 1T05, to form a quaternary complex. After the energy minimization, the geometric quality of the modeled structures was evaluated by PROCHECK32. The protein structures were visualized with the PyMOL molecular viewer33.

Molecular docking

The molecular docking study was performed with the Glide program (version 5.5)34. It should be noted that ADV is a prodrug, and we modeled its bioactive form, ADV diphosphate (ADV-DP), which is generated through stepwise phosphorylation after the oral administration of ADV. The 3D conformation of ADV-DP was constructed according to the coordinates of the ligand, TDV-DP, which was deposited in the template structure of HIV-RT, by deleting a methyl group.

After the ligand preparation, each modeled protein was prepared with the protein preparation wizard in Maestro34. The molecular docking was performed using standard precision (SP) protocols with default parameters. The docking poses were ranked by their glide scores, and the best predicted conformation in each system was used for the subsequent molecular dynamics simulation and binding energy analysis.

Molecular dynamics

Classical molecular dynamics simulations of docked conformation-HBV polymerase complexes were performed using the AMBER 9.0 suite35. To make a reliable comparison, the two systems were treated with the same procedures described below.

The electrostatic potential of ADV-DP was computed with Gaussian03 at the HF/6-31G* level36, then the charges were assigned using the RESP (restrained electrostatic potential) methodology37. The charges and force field parameters for ADV-DP were generated by an antechamber38. The polyphosphate parameters of ADV-DP in our system were identical to those developed by Meagher et al39. The AMBER ff03 force field40 and the general AMBER force field (GAFF)41 were chosen to create the potential of the proteins and ADV-DP, respectively. The systems were then solvated in a truncated octahedral box of TIP3P water molecules 10 Å away from the protein42. Counter-ions of Cl− were then added to obtain the electrostatic neutrality of the systems.

Prior to the MD simulations, stepwise energy minimizations were utilized to relax the two systems by the sander program in AMBER 9.0. First, 3000 minimization of water molecules/ions steps were performed using a combination of steepest descent and conjugate gradient methods. The solute atoms were frozen with a constraint of 100 kcal·mol-1·Å-2. The systems were then relaxed with a constraint of 10 kcal·mol-1·Å-2 for an additional 3000 minimization steps. Next, the side chains of the proteins were minimized with an additional 2500 steps. Finally, the entire systems were minimized with 2000 steps without any constraints.

In the MD simulations, the Particle-Mesh Ewald (PME) method was utilized to address the long-range electrostatic interactions and a cutoff of 10 Å was applied43. All chemical bonds involving hydrogen atoms were constrained using the SHAKE algorithm44, and the time step used for the MD simulations was 2 fs. The system was gradually heated from 0 to 300 K over 50 ps molecular dynamics in the NVT ensemble. Subsequently, a short dynamics of the entire system was performed at 300 K for a further 110 ps to equilibrate the solvent. Finally, 3.7 ns MD simulations were conducted with a periodic boundary condition in the NPT ensemble (300 K, 1 atm). The analyses of MD trajectories were performed with the Ptraj module of AMBER 9.0.

MM-PBSA calculations

The binding free energy (ΔGbind) of ADV-DP in the two systems was calculated using the MM-PBSA method45. ΔGbind was computed as follows:

ΔEgas, the sum of the molecular mechanical energies, is composed of three parts: internal energy (ΔEint), electrostatic potential (ΔEele) and van der Waals potential (ΔEvdw). The solvation free energy (ΔGsolv) can be divided into contributions from the polar solvation free energy (ΔGPB) and the nonpolar (ΔGSA). In the PB calculation, the solvent and solute dielectric constant values were set to 80 and 2, respectively. The LCPO method46 was used to estimate the non-polar contribution as follows: ΔGSA=0.0072·ΔSASA, where a probe radium of 1.4 Å was set for the SASA (solvent-accessible surface area) calculation35. Because of the high computational demand for large molecules, such as proteins, and the rough estimate of the entropic contribution (−TΔS)47,48 and because we are comparing the relative binding affinities of two similar systems (ie, HBV polymerase genotypes B and C complexed with the same inhibitor, ADV-DP), the entropic contribution of both systems is disregarded in this study. The energy difference between the two systems (ΔΔGC–B) is as follows:

This equation was used to measure the sensitivity of genotype C over B with regard to the interaction with ADV.

In our work, 100 snapshots extracted from the last 800 ps simulation trajectories of each system were used to calculate the binding free energy of ADV-DP in the HBV polymerase-DNA-Mg complex.

Free energy decomposition and computational alanine scanning

The contribution of each residue in the HBV polymerase was then analyzed using the MM/GBSA decomposition procedure49,50. The ΔGGB was calculated using the generalized Born (GB) model51, and the ΔGSA was computed with the ICOSA method49. All the energy components were calculated using the same 100 snapshots.

Next, the computational alanine scanning was performed for the two complexed systems. Single site alanine mutations at residues 236 and 238 in both genotypes B and C proteins were generated by deleting atoms and truncating the mutated residues at Cγ and replacing them with a hydrogen atom52. The binding free energies of the alanine mutants were calculated using the MM-PBSA approach described above. The same 100 snapshots obtained from the last 800 ps of the MD trajectory were used.

Results

Homology model of HBV polymerase

A total of six standard reference sequences of genotypes B and C were retrieved using the National Center for Biotechnology Information (NCBI) genotyping tool (Supplementary Information Table S1). The sequence analyses of these reference sequences revealed that amino acid (AA) 238 located near the ligand binding pocket is different among the six consensus genotypic reference sequences (ie, Gln/His/Asn is located at position 238 in genotype B and Asn is located at position 238 in genotype C), and all the other polymorphic sites between the two genotypes were found distant from the binding pocket. Further statistical analyses of the published sequences reported in NCBI confirmed the significant difference at position 238 between genotypes B and C. As shown in Table 1, the occurrence rates of His and Asn at this position are 81% in genotype B and 86% in genotype C, respectively.

Table 1 Statistics of a few major amino acid types at position 238 for different HBV genotypes.
Full size table

Based on the above analyses, the amino acid at position 238 was set to His for genotype B and Asn for C. The overall consensus reference genotype sequences used to construct the homology models of the HBV polymerase genotype B and genotype C share a high identity of 93.3% (Supplementary Information Table S2). The sequence alignment of the HBV polymerase with HIV-1 RT is shown in Figure 2. The final homology models of the HBV polymerase genotypes B and C are shown as ribbon diagrams in Figure 3. Similar to the HIV-1 RT, the modeled structures of the HBV polymerase also have three typical subdomains: finger (1 to 49 and 90 to 172), palm (50 to 89 and 173 to 267) and thumb (268 to 344)25 (Figure 3B).

Figure 2
figure 2

Alignment of HBV polymerase and HIV-1 RT sequences showing the identical (dark shaded) and closely related (light shaded) residues.

Full size image
Figure 3
figure 3

(A) Overall view of superposition of template structure 1T05 (magenta ribbon) with HBV polymerase genotype B (cyan) and C (green) model structures from homology modeling, ligand ADV-DP and two Mg2+ ions are highlighted in magenta and gray. (B) The 3D structural model of finger, palm, and thumb subdomains in HBV polymerase are represented in green, blue and salmon ribbons, respectively. DNA template-primer duplex is shown in orange cartoon mode. (C) The structurally conserved catalytic subdomains are displayed in colored ribbons: Domains A (cyan), B (yellow), C (pink), D (green), E (blue), F (magenta), and G (marine).

Full size image

As illustrated in Figure 3A, the two HBV polymerase models and HIV-1 RT structure resemble to each other in the palm and thumb regions; however, there is a major difference in the finger subdomain. The HBV polymerase has a larger finger region than the HIV RT and contains three long loops of 27, 11, and 9 residues in the HBV finger domain. Overall, the high similarity can be observed between the HIV-1 RT X-ray structure and the modeled HBV polymerase, particularly in the conserved catalytic subdomains, A–G53,54, which are highlighted by different colors in Figure 3C. The structural variations of those two genotypes existed primarily in the flexible loop domain that is distant from ligand active site.

The homology models of the HBV polymerase genotypes B and C were then evaluated by PROCHECK32. The overall G factors were -0.2 and -0.22 for genotypes B and C, respectively, indicating that the molecular geometries of the constructed models are stereochemically acceptable. The Ramachandran plots of the two models are shown in the supporting information (Supplementary Information Figure S1). Of all the residues, 81.7% and 82.2% are located in the most favored regions for genotypes B and C, respectively. Further analysis of the models revealed that all the disallowed residues except Met204 located in the YMDD motif were located distant from the putative ligand-binding site of the HBV polymerase. Consistent with Daga et al25, the corresponding Met residue was also found in the disallowed region of many HIV-1 RT crystal complexes, suggesting that the unusual configuration of Met204 is a valid and conserved feature. These evaluation results show that our homology models of HBV polymerase are of reasonable quality.

Docking studies of ADV diphosphate (ADV-DP) with HBV polymerase genotypes B and C

The superposition of the two docked complexes with the crystal protein, 1T05, suggests that the ADV-DP was docked into the active site with favorable poses, both bearing a close resemblance to the pose of the analogous ligand, TDV-DP, in 1T05. The best possible binding modes of ADV-DP in the HBV polymerase active sites of genotypes B and C are illustrated in Figure 4. Both complexes are stabilized by extensive hydrogen bonding networks. The nitrogen base of ADV-DP displayed two hydrogen bonds with the complementary base pair, dTMP, in the template chain. Moreover, it formed a π–π stacking interaction with the DNA base in the primer chain. The β-phosphate and γ-phosphate of the ligand formed a hydrogen bond with the backbone amide NH of Ala86 and Ala87. The side chains of Arg41 and Lys32 were involved in ionic interactions with the phosphate group of ADV-DP. Arg41, in the systems of genotypes B and C, formed three and two hydrogen bonds with ADV-DP, respectively. The carbonyl groups of Asp83, Asp205 and Val84 together with the three phosphate groups of ADV-DP formed a metal chelating interaction with the two Mg2+ions that were present in the active site. The docking scores of the two genotypes, B and C, are also at the same level, with gscores of -10.25 and -10.07, respectively.

Figure 4
figure 4

Binding modes of: (A) HBV polymerase genotype B and (B) genotype C with ADV-DP. Two Mg2+ ions are shown in gray colored spheres. The ligand ADV-DP is shown in magenta colored ball-and-stick mode. Mg2+ coordination and hydrogen bond interactions are shown with black dotted lines.

Full size image

Assessing the binding mode of the different genotypes by MD simulation

To further validate our docked model, a molecular dynamics simulation was preformed to obtain a more precise ligand-receptor state that is close to the natural behavior of the complex. The conformational stability of these two HBV polymerase-ADV complexes was assessed by performing a 3.8 ns molecular dynamics simulation in an explicit water model under the same conditions. The root mean square deviation (RMSD) values with respect to the initial structures were calculated along the 3.8 ns trajectories and are shown in Figure 5. The backbone RMSD values of the two systems tend to be convergent after 2 ns of simulation, with fluctuations approximately 4.5 Å and 4 Å for the structures of genotypes B and C, respectively. The RMSD change of the conserved catalytic subdomains, which is highlighted in Figure 3, exhibits a similar tendency. Intuitively, the lower conserved subdomain RMSD values of genotype C may suggest that the complex of ADV-DP bound within the HBV polymerase genotype C is in a more stable state.

Figure 5
figure 5

Root mean squared deviation (RMSD) of HBV polymerase homology model during the molecular dynamics simulation with respect to its starting (docked) pose.

Full size image

Figure 6 summarizes distance fluctuations of the corresponding atom pairs that are involved in the hydrogen bond interactions. As shown, some significant hydrogen bond interactions that were observed in the initial docking model (Figure 4) were maintained during the MD simulation. For example, two hydrogen bonds between the 6-amino purine ring of the ADV-DP and the paired base, dTMP, and the hydrogen bond between Arg41 and the α-phosphate of ADV-DP were maintained through the entire simulation of both systems. Moreover, we may also find a newly formed hydrogen bond between the γ-phosphate of ADV-DP and the hydroxyl group of Ser85 during the MD simulations.

Figure 6
figure 6

Distance graphs of some important hydrogen bonding interactions. Two maintained hydrogen bonds: (a) between N1 of ADV-DP and NH of dTMP; (b) and between NH2 of ADV-DP and carbonyl O of dTMP. Panel (c) shows distance variations between α-phosphate of ADV-DP and the residue Arg41. (d) The residue Ser85 was distant in the beginning of the MD simulation; loop movement pushed the side chain amide toward the γ-phosphate of ADV-DP as described. (e) Distances between the hydroxyl group of Ser85 and NH2 of Asn236. (A) Genotype B; (B) Genotype C.

Full size image

The simulation system of genotype C showed superiority during the shorter time needed to form a steady hydrogen bonding interaction. Furthermore, the side-chain conformation of Ser85 in the genotype C system was further stabilized by a hydrogen bond with Asn236 in a loop region (Figure 6E), which was consistent with the binding pattern reported by Langley et al23. However, in the genotype B system, the hydrogen bonding interaction between Ser85 and Asn236 was absent. The distance fluctuations that were observed from the MD simulations suggest that the different AA compositions at position 238 may affect the hydrogen bond formation between Ser85 and Asn236.

MM/PBSA free energy analyses

The binding free energies obtained from the MM/PBSA method are shown in Table 2. The table also lists the detailed contribution of the various energy components. It is clear from the table that the ADV-DP showed better binding affinity towards HBV polymerase genotype C than genotype B (ΔΔGC–B=-20.96 kcal/mol).

Table 2 MM-PBSA binding free energy estimations (kcal/mol) of HBV polymerase-DNA-Mg2+-ADV-DP complex of genotype B and C in the MD simulations.
Full size table

To provide quantitative information regarding the detailed contribution of the residues in the different HBV polymerase genotypes, free energy decomposition analysis was utilized to decompose the contributions to the binding energies on a per-residue basis. Many residues were found to exhibit influences on the binding to the inhibitor, ADV-DP. Several amino acid residues, Lys32, Arg41, Asp83, Ser85, Ala86, Ala87, Phe88, Asp205, and Asn236, and two DNA base pairs were validated to have significant positive contributions with an absolute energy greater than 1 kcal/mol. This is highlighted in Figure 7. Notably, ADV-DP binds to similar key residues in the two systems, which is consistent with the results of the docking study and MD simulation. The favorable contributions of Lys32, Arg41, Asp83, Ser85, Ala86, Ala87, and dTMP are due to the hydrogen bonding interactions with ADV-DP. The purine ring of ADV-DP formed a strong π-π stacking interaction with the primer DNA base, which also showed an obvious hydrophobic interaction with Phe88. Phe88 is located within the hydrophobic pocket comprised of Ala87, Phe88, Ile180, and Met204. This finding is consistent with the report by Daga et al25. In addition, for most of the key residues, a slightly stronger binding energy contribution was found in the genotype C system.

Figure 7
figure 7

Decomposition of binding energy (GBTOT) on a per-residue basis, only residues making significant favorable or unfavorable contribution were shown (|GBTOT|≥1.0 kcal/mol).

Full size image

As observed in Figure 7, positive contributions were made for all the key residues except Asp205, which possibly arises from the repulsive interaction between the negatively charged side chain of Asp205 and the negatively charged phosphate group of ADV-DP. By contrast, the Mg2+ions between Asp205 and the phosphate group of ADV in these two stable complexes (Figure 4) mediate a favorable coordination linkage between D205 and ADV-DP. The genotype-dependent, polymorphic residue 238 shows less contribution to the binding energy (≤1.0 kcal/mol). A lower residue contribution to genotype B than C was found at positions 236 and 238, which was verified by the computational alanine scanning results (shown in Table 3). Both the alanine mutations at residues 236 and 238 lead to losses of sensitivity, as reflected by an increase of the binding free energy difference between genotypes C and B (ΔΔGC–B). For the N236A mutant, the ΔΔGC–B increased from -20.96 (ΔΔGC–B, wild type) to −18.66 kcal/mol (ΔΔGC–B, N236A). For the H/N238A mutant, the ΔΔGC–B increased to -19.96 kcal/mol (ΔΔGC–B, H/N238A). These results confirmed our hypothesis that residues 236 and 238 contribute to the sensitivity of genotype C over B regarding the interaction with ADV-DP. Moreover, Asn236 is more important because the alanine mutation at this residue site causes a more significant loss of sensitivity.

Table 3 Changes of MM-PBSA binding free-energy (kcal/mol) caused by the alanine mutations at the residue 236 and 238 in two genotypes systems.
Full size table

Discussion

Hepatitis B is a growing problem worldwide. Genotypes B and C are the two predominant types of HBV in Asia. Clinical studies have been conducted to investigate the influence of specific genotypes on ADV treatment; however, the results are not yet conclusive. To further explore the underlying nature of the sensitivity of different genotypes of HBV polymerase to ADV, we used a computational protocol that combined homology modeling and binding free-energy calculations with MD simulations.

For the homology modeling of different genotypes of HBV polymerase, the representative sequence selection is of great significance. The statistical sequence analyses shown in Table 1 reveal that residue 238 near the binding pocket is not only a polymorphic site but is also a genotype-specific site (ie, His238 is present in genotype B and Asn238 is present in genotype C). This observation is consistent with previously reported sequence analysis results. An investigation of 192 treatment-naïve Chinese patients showed that position 238 is a genotype-dependent polymorphic site, in which 52/55 genotype B patients have a His at this position and 131/137 of genotype C patients have an Asn55. Moreover, Rhee et al demonstrated that the AA types His238 and Asn238 are typical for B and C genotypes, respectively9.

In the polymerase/RT region, the sequence identity among various HBV strains is significant; however, the identity between HBV and other viruses is relatively low. HIV-1 RT shares approximately 25% sequence identity with HBV22 and has available crystal structures and is one of the nearest relatives of HBV polymerase. Many of the functionally important residues that are involved in the key protein-ligand interactions in HIV-1 RT are conserved between the HIV-1 RT structure and the HBV polymerase. These include the YMDD motif, the MGY motif and the catalytic triad of aspartic acid residues. For example, the catalytic triad Asp83, Asp205, and Asp206 of HBV polymerase correspond to Asp110, Asp185, and Asp186 in HIV RT (Figure 2). These three Asp residues are responsible for the coordination of two Mg2+ ions, which are very important for the activity of the polymerase. These data suggest that HIV-1 RT shares sufficient structural and functional similarity with HBV polymerase to serve as an excellent template to build a 3D structure for the HBV polymerase22,23,24. In our work, the sequence alignment of our homology model referred to the alignment of the latest sequence published by Daga et al. Based on this sequence, the resulting homology model can match the results of the experimental mutation studies of the HBV polymerase25. Only a minor modification was made at site Asn236 according to the work of Das et al22, in which Asn236 of the HBV polymerase was aligned with Pro217 of HIV-RT instead of the Thr216 as in Daga's work. Asn236 is an important residue in the binding pocket of HBV polymerase and it is a primary drug resistance mutation site to ADV-DP56. However, the alignment of Daga et al would result in the Asn236 in the constructed structure to be located too distant (∼7 Å) from the nearest atom of the ligand (CTP), as suggested by their docking study. Therefore, for this residue, we chose the alignment of Das et al and used residue Pro217. As a result, the position of Asn236 in the constructed homology model is more feasible and is located in close vicinity to the binding site of ADV-DP. Apart from this, there are no other changes.

After the three-dimensional structure was successfully constructed, the ADV-DP was docked into the HBV polymerase-DNA complex. From the binding models, similar hydrogen-bond networks are observed for the ADV-DP in the HBV polymerase genotypes B and C. In addition, the gscores of the two systems are at the same level. All these indicate that the binding affinity that was predicted by docking could not discriminate between genotypes B and C. We must also note that molecular docking cannot fully account for the impact of protein flexibility, and a small change in receptor conformation may lead to a significant variation of binding affinity. Because identical structure templates were used in the homology modeling of the HBV polymerase genotypes B and C, the similar binding patterns of the two systems may simply be a result of the similar static protein conformations that were modeled. Therefore, to consider the flexibility of both the ligand, the protein and also their interactions with the surrounding water molecules, molecular dynamics simulations in combination with MM/PBSA free energy calculations were performed.

In summary, the calculated results support the hypothesis that the HBV polymerase genotype C is more sensitive to ADV than genotype B. Moreover, residue N236 and the polymorphic site 238 were further investigated and shown to play an important role in contributing to the sensitivity of genotype C over B in the interaction with ADV-DP. We believe that this work should provide a basis for the further investigation and understanding of the effect of polymerase genotypes on ADV treatment. This allows for the optimization of individual treatment and the minimization of side effects associated with HBV therapy.

Author contribution

Ming-yue ZHENG and Ai-long HUANG designed the study, supervised the research and revised the manuscript. Jing LI performed the research, analyzed the data and wrote the manuscript. Yun DU, Xian LIU, Qian-cheng SHEN, Xiao-min LUO, and Hua-liang JIANG discussed the results and commented on the manuscript.