Molecular docking and dynamic simulation analysis of Hepatitis E virus protease in complexing with the E64 inhibitor

Abstract The unavailability of a suitable treatment for human Hepatitis E virus (HEV) infection necessitate the development of anti HEV drugs. The HEV papain-like cysteine proteases (HEV PCP) is a crucial target to prevent viral replication and progression. E64 is a known HEV PCP inhibitor; however, its molecular mechanism of inhibition is not yet known. Since the crystal structure of HEV PCP is not available, the primary focuses of the present study was to refine the predicted HEV PCP structural model by molecular dynamics (MD) simulation. Further, we performed a 200 ns MD simulation to understand the structural complexity of HEV PCP and the effect of E64 binding with HEV PCP. The E64 binding with active site residues Gln48, Thr51, Gln55, Cys52, Ser81, Gln 98, Cys 132, Arg158, His159, Asn 160 and Ala96 leads to reduced fluctuations in the residue at N-terminal (18–41) that include the CHC motif (26–28). However, most of the other non interacting residues, including the inter-domain linker region (46–87), showed increased fluctuations in the HEV PCP-E64 complex. The residue Asp21 and Ala96 are involved in the formation of interdomain interactions in the HEV PCP apo enzyme. While in the PCP-E64 complex, E64 binds to Ala96 and creates a steric hindrance to prevent interdomain interactions. Thus, the E64 binding reduces interdomain interactions and restrict domain movements in the HEV PCP-E64 complex. This information will be important for the chemically designing more effective derivatives of E64 developing HEV PCP specific inhibitors. Communicated by Ramaswamy H. Sarma


Introduction
HEV mediated Hepatitis E liver disease caused approximately 44,000 deaths in 2015, accounting for 3.3% of the mortality due to viral hepatitis (WHO, 2021). In addition, a very high mortality rate has been observed in HEV infected pregnant women in their third trimester (Wu et al., 2020). Although human HEV infection has been found worldwide, the disease is most common in East and South Asia. Currently, chronic HEV infection, HEV-related acute hepatic failure, and extrahepatic manifestations caused by HEV infection have been frequently reported (Himmelsbach et al., 2018;Ju et al., 2020;Mikulska et al. 2021;Nan et al., 2017).
Current therapeutics used to treat HEV infection, including antiviral agents ribavirin and Interferon-a, have severe side effects, are contraindicated in pregnant women, reported treatment failure (Kenney & Meng, 2015). More than 20 second-generation direct-acting antiviral agents (DAAs) are in phase II/III clinical trials. However, their efficacy against HEV is doubtful as most of them have been specifically designed against hepatitis C viral protein that does not share high sequence similarity with hepatitis E virus (Kinast et al., 2019). Currently, an HEV vaccine has been developed and is licensed only in China, which would not help treat immunocompromised individuals (Kenney & Meng, 2015).Thus, there is a need to develop affordable, potent, novel, safer and more effective antivirals to treat severe hepatitis E diseases in high-risk populations.
HEV, a 7.2-kb quasi-enveloped, positive-sense, singlestranded RNA virus, belongs to the family Hepeviridae (Wang & Meng, 2021), which contains three partially overlapping reading frames ORF1 to ORF3. All non-structural proteins (nsPs) required for viral replication have encoded by open reading frame1 and are transcribed as a single transcript. ORF1 encodes methyltransferase, RNA-dependent RNA polymerase and RNA helicase proteins. ORF1 also contains a few poorly understood domains, including a hypervariable region, a putative protease, X and Y domains (LeDesma et al., 2019) Recently putative protease domain has been expressed in Ecoli and baculoviruses. Mutation in conserved cysteine and histidine residues in the putative protease inhibits proteolytic activity showed its role in ORF1 processing and viral replication (Paliwal et al., 2014;Parvez, 2013;Parvez & Khan, 2014).
Most RNA viral proteases are categorized either as chymotrypsin-like or papain-like proteases (PCP). Chymotrypsin like proteases, presents three catalytic triads, Cys-His-Asp, Cys-His-Glu and Ser-His-Asp residues. In PCPs, a conserved Cys-His dyad is assisted by unique Asn and Asp residues. In both the chymotrypsin-like or PCPs, Zn 2þ is located at the side opposite the active centre. The PCPs have been characterized by the presence of two domains, an a-helix and a b-sheet domain forming a deep cleft containing Cys-His-Asn residues along with Gln acting as a substrate-binding site (Saraswat et al., 2020).
The HEV-protease has been reported acting as PCP with conserved catalytic dyad (Cys and His, but it has also been reported as a chymotrypsin-like protease (Paliwal et al., 2014;Parvez, 2017). In our previous study, the predicted 3 D model showed papain-like two domain architecture (N-terminal helical domain and a C-terminal b-sheet domain) Cys52-His159-Asn160 residues forms the active site triad. The presence of residue Gln48 and Cys-His and Asn triad at the active site brings HEV-protease closer to PCPs. The predicted active site for substrate binding was present between the two domains.
Structure-based drug design has emerged as a powerful method for designing novel potent inhibitors for the target enzymes with known 3 D molecular structure (Lu et al., 2018). Since the 3 D crustal structure of HEV PCP is not available, the present work focused on refinement of the predicted HEV PCP structural model (100 ns) and further studied the effect of HEV PCP-E64 docking by MD simulations (200 ns).
Ligand restriction grid box was set (center ¼ (41.207, 43.627, 44.299 Å), size ¼ (40, 40,40 Å)) covering the entire structure. The exhaustiveness set 50; this allows for a global search of the ligand E64 to dock to a more favourable place on the protein. Finally, nine different conformations were generated. The binding interactions between the HEV PCP structure with ligand E64 were plotted using ligand interaction diagrams of Maestro (Maestro, 2017) and LigPlot þ v.2.2 (Laskowski & Swindells, 2011).
The topology of HEV PCP was generated using pdb2gmx modules of GROMACS. First, the predicted model was placed in a cubic box using the editconf module, keeping protein HEV PCP at the centre and 10 Å distance from the edges. Next, the system was solvated with a simple point charge water model (SPC216) to attain real dynamics. Finally, solvated systems were energy minimized using the Steepest Descent algorithms with 5000 steps cutoff. Further, individual systems were slowly warmed up using a V-rescale thermostat with a coupling constant to reach 300 K to perform equilibration in NVT (Number of atoms, Volume of the system, and Temperature of the system) ensemble. Later, the solvent density was maintained using a Parrinello-Rahman barostat with a coupling constant at 1 bar and 300 K to perform equilibration in NPT (Number of atoms in the system, the pressure of the system and temperature of the system) ensemble.
With restrained position dynamics, the systems were equilibrated for about 100 ps under both ensemble processes (NVT and NPT). The LINCS algorithm has been used for constraining all the bonds. Finally, the systems were submitted to MD simulation for 200 ns to observe the stability of HEV PCP. The MD trajectories were saved every 10 ps with a time step of 2.0 fs. GROMACS and visual molecular dynamics (VMD 1.9.1) (Humphrey et al., 1996) have been used to analyze the resultant trajectories. All the graphs have been plotted using Grace 5.1.23 program (Turner, n.d.). One more MD simulations were performed with Desmond, as implemented in the Schr€ odinger package (Shaw, 2021).

Structural refinement
The 3D structure of HEV PCP developed by homology modelling has been refined by 100 ns MD simulations. During MD simulations, no major change in the energy was observed, and RMSD changes have also stabilized within the first five ns of simulations ( Figure S1). The Ramachandran plot analysis shows that the structure model stereochemical quality increased considerably by MD simulation ( Figure S2). After MD simulation, the residue count increased significantly in the favoured region (82.6%-90.6%); decreased in the outlier (1.4%-0%) and permitted regions (19.9%-9.4%) ( Figure S2(B) and Table S1).
Furthermore, validation of the refined structure model, based on Z-score (for overall model quality) and energy of residues (for local model quality), indicates that the model is of good quality and are energetically stable. The Z-score of the refined model was À3.95 (in the range of Z-score for PDB structures determined by NMR) ( Figure  S2(C)). Analysis of the 3 D model based on residue energies indicates that most residues are in a low energy state ( Figure S2(D)). Energy versus sequence length plots for a window size of 10 and 40 amino acids further validate the refined model's stability. Positive energies for the N-terminal sequences were observed for the window size of 10 and 40 amino acids, whereas no positive energies for C-terminal sequences were observed ( Figure S2(E)).
Potential energy versus time plot was drawn for the trajectories captured at an interval of 10 ps. Only minor energy fluctuations were observed during the MD simulation with minute energy change, indicating no major changes in the protein conformation ( Figure S3(A)). Secondary structure analysis of initial and refined HEV PCP structures showed no major change in the protein structure ( Figure S4).

Molecular docking of E64 with HEV PCP
The HEV PCP structure cavity showed Cys50, Cys52, His159, Asn68, Ser81, Asn154, Asn160, and Gln48, Gln55, Gln98 residues line the predicted active site. The Catalytic tried (Cys, His and Asn) is contributed by two domains and a highly conserved glutamine residue (Gln48). E64, the known HEV PCP inhibitor, was docked with refined HEV PCP structure using blind docking (Figure 1, Table S2). The observed binding energy of E64 with the HEV PCP complex was À6.4 (kcal/ mol). Usually, a compound is predicted to have activity against a protein when binding energy is less than À6.0 (kcal/mol). The E64 binding site contains residues of catalytic dyads such as Cys 52 and His159 (Table S2). Further, the Ligand interaction diagrams of the HEV PCP-E64 complex showed binding within the active binding site residues Asp21 and Ala96 (Figure 1, Table S2).

RMSD analysis
RMSD values varied between 0.43, 0.54, and 0.56 nm for HEV PCP and HEV PCP-E64 complex and replica of HEV PCP-E64 complex, respectively. An overlapping plot between HEV PCP and HEV PCP-E64 complex showed minor RMSD value changes (Figure 2(A)). In both HEV PCP and HEV PCP-E64 complex, observed that the system remained steady till 200 ns with minor fluctuations with time throughout the simulation. The HEV PCP-E64 complex showed more fluctuations as compared to that of the HEV PCP structure. A little higher RMSD in HEV PCP-E64 complex form was observed after a 68 ns period and finally stabilized at 100 ns. RMSD values of both the structures, the HEV PCP, were more stable than that in the HEV PCP-E64 complex, suggesting that inhibitor induces minor structural changes in HEV PCP structure.

SASA analysis
A comparison of the analysis of global and local stability and compactness of HEV PCP and HEV PCP-E64 complex during simulation has been summarized in Figure 2(B). SASA was investigated for both the systems and observed a similar area has exposed to the solvent except for a brief period (65-120 ns) where comparatively increases in SASA was observed during the MD simulations. The HEV PCP-E64 complex showed a relatively similar area exposed to the solvent after 120 ns of MD simulations (Figure 2(B)). Average SASA values varied between 106 and (107, 111) nm2 for HEV PCP and HEV PCP-E64 complex, respectively.

RMSF analysis
All structures simulations have been used for the RMSF analysis (Figure 2(C)). The peaks RMSF plot indicates regions of the protein that fluctuated the most during the simulation. In addition, it has been observed that HEV PCP-E64 complex residues fluctuations were higher than HEV PCP, which clearly shows that the backbone atoms undergo more fluctuations upon inhibitor binding.
The RMSF plot of the residue's backbone indicates a similar pattern of fluctuations for both MD simulations ( Figure  2(C)). Minor fluctuations have been observed in the terminal, turn, and loop regions; whereas, major fluctuations have been seen in the loop (residue no. 26-28, respectively). The largest fluctuation of 0.6 nm was observed for His 27 residue in the largest loop inter-domain linker region followed by (Residue no. 46-87) that connects the a helix 1 to the b1 strand ( Figure S4). The third major fluctuation was observed in the loop region of HEV PCP-E64 complex residue Asp 127 connects the b2 strand to the b3 strand. It indicates that E64 binding with active site atoms reduces fluctuations in the residue at N-terminal (18-41), including the CHC motif (26-28). However, most of the remaining residues, including the inter-domain linker region (46-87), showed increased fluctuations in HEV PCP-E64 complex.

Rg analysis
Furthermore, complex compactness was identified by analyzing Rg (backbone vs time). The HEV PCP-E64 complex shows a slightly higher distance than the HEV PCP (Figure 3(A,B)). The Rg plot indicates that the degree of compactness was low initially (high radius of gyration) but changed due to the ligand's slight movement in the pocket of HEV PCP. The system's overall Rg has been observed to be higher than the HEV PCP (Figure 3(A,B)).

H-bonds analysis
Variations in the total number of H-bonds in HEV PCP and HEV PCP-E64 complex during MD simulation have been presented in Figure 3(C,D). The average H-bond per time frame observed was 87, (92, 84), and (2.5,4) for HEV PCP, HEV PCP-E64 complex and inhibitor E64, respectively (Figure 3(C,D)). The number of H-bonds versus time (ns) plot reveals a

PCA analysis
Atomic-level MD simulation offers to explore the overall structural movement in HEV PCP and HEV PCP-E64 complex. Therefore, we performed PCA to identify large-scale collective motions of the HEV PCP and HEV PCP-E64 complex trajectories generated by 200 ns MD simulations (Figure 4).
PCA is often used to characterize different protein conformational dynamics involved in protein folding, open-close mechanism of ion channels. The eigenvalues were derived for both the systems and the first ten hits were considered for analysis. The data showed 84.5%, 88.9% and 91% motions for HEV PCP, HEV PCP-E64 complex and replica of HEV PCP-E64 complex, respectively (Figure 4). The active site residues are in black. The interconnecting loop region starts from Leu46 to Ser87 and is shown in green color. The connecting residues are responsible for movement represented by the dotted black line and connecting residues (Asp21 and Ala94) are shown in blue color.
The figure shows that HEV PCP and HEV PCP-E64 complex exhibited the motions. Then, the first two eigenvectors were selected and plotted against each other (PC1 vs PC2). In this plot, that HEV PCP showed multiple clusters (Figure 4).
The structure of HEV PCP undergoes conformational changes due to the structure constituted by the N and C terminal domains connected by a long loop residue number 46-87 (a stretch of 41 amino acid residues), and there are no disulfide bridges. Hence, further HEV PCP and HEV PCP-E64 complex analyses have been made by focusing on the active site region. Snapshots have been taken after PCA of the HEV PCP, showing transition through the MD simulation ( Figure 5).
Multiple open, semi-open and closed conformations have been observed; these confirmations have been visited numerous times. On the other hand, these confirmations show anticorrelated motions (negative-positive directions). There is a considerable concerted movement in the N-terminal segment of both HEV PCP and HEV PCP-E64 complex (Figure5). Close confirmations were absent in the HEV PCP-E64 complex. It also shows that there is a possibility that the conformational changes (open to close and vice-a-versa) occurred in HEV PCP due to the anticorrelation movement between the PCP and methyltransferase domain, and it did not affect the active site as reflected in Figure 5, Table S3 and Figure S7.
For ligand interaction diagrams analysis, twenty-six snapshots are taken for every eight ns time intervals over 200 ns MD simulation trajectory ( Figures S8-S11). The active site residues have been highlighted in the red colour circle, and residue involved in the loop movement in blue.

Conclusion
HEV PCP plays an essential role in viral replication and propagation; thus, inhibiting this protein would be therapeutic. Structural biologists are continuously resorting to MD simulations to identify a better drug candidate and other therapeutic agents. Here, we refined HEV PCP structural properties using 100 ns MD simulation. The Rg and SASA analysis exposed higher and lower fluctuations, suggesting that protein might have gone folding state; the degree of fluctuation in Rg indicates less rigidity or compactness due to protein folding or ligand binding. In addition, RMSF analysis could bring out some important information on atomic or residual flexibility that might contribute to different conformational states.
The PCA analysis results of HEV PCP and HEV PCP-E64 complex show conformational changes and adopt an alternate conformation like open, semi-open, and closed confirmations. Interestingly, the active site residue interactions remained stable in HEV PCP, even though it undergoes large conformational changes. PCA analysis and conformational clustering supported the above findings. Together, the results are helpful to confirm the conformational change that undergoes in HEV PCP.
This study concludes with a very primary understanding of HEV PCP protein. Although a few significant deviations were marked from RMSD, those did not remain afterwards. This study found a significant change in residual flexibility and its contribution to different possible conformers, as explained in RMSF and PCA analysis. This study may be used in structure-based drug discovery against HEV PCP. In this regard, a study provides a clue for designing novel drugs with better specificity and affinity. As per our MD simulation and docking analysis, the active site residues Gln 48, Thr 51, Gln55, Cys52, Ser81, Gln 98, Cys 132, Arg158, His159, Asn 160 and along with Ala96, which is involved in the loop movements, demonstrated crucial interactions with HEV PCP. The work is significant as it emphasises the modelling of HEV PCP structure and virtual screening of drugs that may be used as scaffolds to develop more potent HEV inhibitors.