HSV-1 Glycoprotein D and Its Surface Receptors: Evaluation of Protein–Protein Interaction and Targeting by Triazole-Based Compounds through In Silico Approaches

Protein–protein interactions (PPI) represent attractive targets for drug design. Thus, aiming at a deeper insight into the HSV-1 envelope glycoprotein D (gD), protein–protein docking and dynamic simulations of gD-HVEM and gD-Nectin-1 complexes were performed. The most stable complexes and the pivotal key residues useful for gD to anchor human receptors were identified and used as starting points for a structure-based virtual screening on a library of both synthetic and designed 1,2,3-triazole-based compounds. Their binding properties versus gD interface with HVEM and Nectin-1 along with their structure-activity relationships (SARs) were evaluated. Four [1,2,3]triazolo[4,5-b]pyridines were identified as potential HSV-1 gD inhibitors, for their good theoretical affinity towards all conformations of HSV-1 gD. Overall, this study suggests promising basis for the design of new antiviral agents targeting gD as a valuable strategy to prevent viral attachment and penetration into the host cell.


Introduction
Herpesviridae is a large family of enveloped double-stranded DNA viruses that includes eight human herpesviruses (HHV) divided into three subfamilies: alphaherpesvirinae (herpes simplex virus [HSV] and varicella-zoster virus [VZV]); betaherpesvirinae (cytomegalovirus [CMV], HHV-6, and HHV-7); and gammaherpesvirinae (Epstein-Barr virus [EBV] and HHV-8) [1]. All herpesviruses, especially those belonging to the alphaherpesvirinae subfamily, are able to establish latent infections that can be reactivated by endogenous or exogenous stimuli, causing clinical symptoms [2]. HSV infections are widely spread around the world, thus representing a considerable public health issue. There are two main types of HSV that infect humans, HSV-1 and HSV-2. Both are responsible for genital herpes, although HSV-1 mostly causes oral infections and is acquired during childhood [3]. HSV-1 infection can also lead to complications, such as keratitis, encephalitis, meningitis, and systemic disease in neonates and immune-compromised patients [4,5]. Currently, antiviral drugs available for the treatment of herpesvirus infections include acyclovir (ACV), penciclovir, valacyclovir, famciclovir, foscarnet, and cidofovir, which terminate viral DNA synthesis by inhibiting the viral DNA polymerase [6]. However, these drugs are not effective in the complete eradication of the infection, but only in reducing the frequency and duration of the episodes. In addition, both the emergence of drug-resistant 1.
Herpesvirus entry mediator (HVEM), a member of the tumor necrosis factor receptor superfamily (TNFR) expressed on activate lymphocytes and in other human tissues including kidney, lung, and liver; 2.
Due to the crucial role that gD plays in host cell fusion, a thorough knowledge of its structural features and of the molecular mechanisms regulating its activity may be useful to develop new inhibitors.
The current study aims at a deeper insight into gD and its cellular receptor interactions by means of computational methods, taking advantage of the HSV-1 gD experimentally solved structures in complex to HVEM and Nectin-1. The complexes resulting from protein-protein docking experiments were submitted to molecular dynamic simulations (MDs). GBPM and free energy calculation analyses revealed the most stable complexes and their pivotal points of interaction. Considering that protein-protein interactions (PPI) represent attractive targets for drug design, a structure-based virtual screening (SBVS) was performed taking into account the key residues useful for gD to anchor human receptors.

Results and Discussion
In order to investigate the binding affinity of gD for its surface receptors and to gain useful structural insights regarding the pivotal interactions occurring at the gD interface, we performed a knowledge-based protein-protein docking employing HADDOCK 2.4 web-server [36]. Due to the involvement of gD both in cell adhesion function and viral entry mechanism depending on its cellular receptors, we considered the X-ray structure of gD both in complex to HVEM (PDB code: 1JMA) and Nectin-1 (PDB code: 3U82), respectively at 2.65 and 3.16 Å resolution. The complexes were prepared using Protein Preparation Wizard and uploaded in PDB format to the HADDOCK server. The input docking parameters were discussed in the Section 3.1. All the generated complexes were clustered according to their HADDOCK score, calculated as a weighted sum of a variety of energy terms (including van der Waals, electrostatic, desolvation, and restraint violation energies) and buried surface area (BSA); the Z-score value was also calculated to select the best cluster with respect to all obtained clusters: the most negative Z-score is indicative of the top ranked cluster. For the gD-HVEM and gD-Nectin-1 complexes, we obtained 188 structures gathered in 9 clusters ( Figure S1), and 198 structures in 2 clusters ( Figure S2), respectively. For both complexes, the docking results of each cluster were reported in Table 2. The best clusters (the lowest in energy) were analysed in terms of interactions using Maestro graphical user interface ( Figure 2) and were aligned to the experimental structures through the Protein Structure Alignment tool.  The best clusters (the lowest in energy) were analysed in terms of interactions using Maestro graphical user interface ( Figure 2) and were aligned to the experimental structures through the Protein Structure Alignment tool. By overlapping ( Figure S3), it resulted that for gD-HVEM and gD-Nectin-1 complexes, the RMSD value between the best docked pose and the X-ray structure was 1.54 and 0.80 Å, respectively. All the residues involved in the PPIs are summarized in Table 3. By overlapping ( Figure S3), it resulted that for gD-HVEM and gD-Nectin-1 complexes, the RMSD value between the best docked pose and the X-ray structure was 1.54 and 0.80 Å, respectively. All the residues involved in the PPIs are summarized in Table 3. Table 3. Residues involved in the interaction of gD with HVEM and Nectin-1.

gD-HVEM
gD-Nectin-1 gD HVEM gD Nectin-1 The best obtained docked structures for each complex were refined using Protein Preparation Wizard and energy minimized with OPLS_2005 force field [37]. In order to explore any potential conformational changes of the gD interface, we submitted 100 ns of MDs for both complexes. The stability of MDs trajectories was monitored by the RMSD trend of the protein's backbone atoms from its initial conformation. The average RMSD values of gD-HVEM and gD-Nectin-1 complexes were 3.05 Å and 2.50 Å, respectively. By monitoring the distances between the interface's residues for gD-HVEM, we observed that several interactions previously detailed in Table 3 were maintained during MDs run, except for Asn15, Gly19, and Lys122 of gD, which were subjected to greater fluctuations, thus preventing the contacts with HVEM residues. Instead, Ala7 gained contacts with Ser20 during the whole MDs ( Figure 3).  Table 3. Residues involved in the interaction of gD with HVEM and Nectin-1.

gD-HVEM gD-Nectin-1 gD HVEM gD Nectin-1 Hydrogen bonds
Lys10 The best obtained docked structures for each complex were refined using Protein Preparation Wizard and energy minimized with OPLS_2005 force field [37]. In order to explore any potential conformational changes of the gD interface, we submitted 100 ns of MDs for both complexes. The stability of MDs trajectories was monitored by the RMSD trend of the protein's backbone atoms from its initial conformation. The average RMSD values of gD-HVEM and gD-Nectin-1 complexes were 3.05 Å and 2.50 Å, respectively. By monitoring the distances between the interface's residues for gD-HVEM, we observed that several interactions previously detailed in Table 3 were maintained during MDs run, except for Asn15, Gly19, and Lys122 of gD, which were subjected to greater fluctuations, thus preventing the contacts with HVEM residues. Instead, Ala7 gained contacts with Ser20 during the whole MDs ( Figure 3).  Regarding gD in complex to Nectin-1, the interactions engaged by Gln27 and Gln132 of gD at the interface with Nectin-1 were lost, but it was found that Tyr38 was located in proximity to Gly86 of Nectin-1 during MDs (Figure 4). Regarding gD in complex to Nectin-1, the interactions engaged by Gln27 and Gln132 of gD at the interface with Nectin-1 were lost, but it was found that Tyr38 was located in proximity to Gly86 of Nectin-1 during MDs (Figure 4). To better characterize the gD interface, a deeper analysis using GBPM was performed on all frames of MDs for each system. This method helps to map the key hotspots responsible for PPI by combining GRID molecular interaction fields (MIFs) according to the GRAB tool algorithm [38]. We considered gD as guest and HVEM and Nectin-1 as hosts. Three GRID probes, such as DRY, N1, and O, were chosen to mimic the hydrophobic, H-bond donor, and acceptor areas, respectively. Taking into account an energy threshold above the 30% from the global energy minimum GRID points, we summarized the pivotal residues up to 3 Å from GBPM points. The contribution of each residue was derived by the summa of its GBPM points energy in the matching frames (see Tables S1 and S2). After calculating the average score based on the total number of frames, the key hotspots were split into quartiles to increase the clarity for the reader: quartile 1 (Q1) includes the residues with the major contribution to PPI until quartile 4 (Q4), which contains residues with the weakest interactions, during the entire trajectories (Table 4). To better characterize the gD interface, a deeper analysis using GBPM was performed on all frames of MDs for each system. This method helps to map the key hotspots responsible for PPI by combining GRID molecular interaction fields (MIFs) according to the GRAB tool algorithm [38]. We considered gD as guest and HVEM and Nectin-1 as hosts. Three GRID probes, such as DRY, N1, and O, were chosen to mimic the hydrophobic, H-bond donor, and acceptor areas, respectively. Taking into account an energy threshold above the 30% from the global energy minimum GRID points, we summarized the pivotal residues up to 3 Å from GBPM points. The contribution of each residue was derived by the summa of its GBPM points energy in the matching frames (see Tables S1 and S2). After calculating the average score based on the total number of frames, the key hotspots were split into quartiles to increase the clarity for the reader: quartile 1 (Q1) includes the residues with the major contribution to PPI until quartile 4 (Q4), which contains residues with the weakest interactions, during the entire trajectories (Table 4).
For both protein-protein complexes, 100 frames extracted by MDs were adopted to calculate the relative binding free energy (∆Gbind) using Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) methodology, as applied in a recent study [39]. The results of the calculated ∆Gbind trend for gD-HVEM and gD-Nectin-1 are depicted in Figure 5. After investigating the PPIs' structural details, we focused on the key residues of gD liable for the interaction with the analyzed human receptors to find potential ligands able to prevent the connection with them. According to the literature [40], several gD N-terminus residues at position 7-15 and 26-29 are responsible for the bind with HVEM, After investigating the PPIs' structural details, we focused on the key residues of gD liable for the interaction with the analyzed human receptors to find potential ligands able to prevent the connection with them. According to the literature [40], several gD N-terminus residues at position 7-15 and 26-29 are responsible for the bind with HVEM, whereas two residues, such as Asp26 and Tyr38, also belonging to gD N-terminus, are relevant for the connection with Nectin-1 [41]. Apart from the shared Asp26 and Gln27 residues, both human receptors interact with different portions of gD. It was underlined that HVEM and Nectin-1 show non-reciprocal competition for binding to gD [42]. Indeed, the interaction between HVEM and gD induces a conformational change in gD that results in the formation of the N-terminal hairpin structure that masks the binding site of Nectin-1 [42]. Similarly, the interaction with Nectin-1 induces a new conformation of the gD N-terminus that prevents the HVEM binding [4]. Considering therefore that gD is not able to bind both receptors at the same time [41], we evaluated the possible conformational changes of gD during MDs. As previously reported [43], MDs trajectories were clustered based on RMSD matrix using the average hierarchical clustering linkage method, obtaining three representative structures of gD. Taking into account that the pivotal interactions of gD with HVEM and Nectin-1 involve distinct residues, we used all three representative structures generated from each complex ( Figure S4)  We calculated the average G-score value for each cluster, aiming to focus on the compounds able to better recognize both Pocket 1 and 2 (Table S3). From the docking results, it clearly emerged that the tricyclic [1,2,3] triazolo [4,5-h] [1,6] naphthyridine moiety showed a lower ability to recognize both the binding pockets compared to the bicyclic triazolo [4,5-b] pyridine derivatives. Furthermore, regardless of the electronic nature of the substituents at the N-3 phenyl ring of the triazolo-naphthyridine core, no significant differences could be observed.
Overall, among the 85 investigated compounds, we focused on the binding mode of the compounds able to recognize all the representative structures of both pockets with average G-score values lower than −5.00 kcal/mol. Accordingly, four triazole-pyridines derivatives, 14 ( Figure S5    In gD-Pocket 2, 14, 16, 83, and 85 engaged several H-bonds with Leu28, Asp30, Asn227, and a π-π stacking interaction with Phe223, probably due to the absence of the N-terminal extension (Figure 7). average G-score values lower than −5.00 kcal/mol. Accordingly, four triazole-pyridines derivatives, 14 ( Figure S5), 16 (Figure S6), 83 (Figure S7), and 85 ( Figure S8) showed a favored energetic profile in complex with all gD conformations. In detail, considering the most populated cluster resulted for gD-Pocket 1, the hydroxyl group of 14 ( Figure 6A), 16 ( Figure 6B), 83 (Figure 6C), and 85 ( Figure 6D) was anchored to Leu25. Moreover, 16 and 83 formed an additional H-bond between the amino group and Leu25.   Desmond [44]. The results of MDs were investigated in terms of stability and conformational flexibility in the presence of the selected compounds. The stability of the complexes was evaluated by calculating the RMSD of the protein's backbone atoms from its initial to final conformation over the whole simulation. By RMSD analysis, we observed that the most promising compounds maintained overall stability throughout MDs in both Pockets, as shown in Figure 8. In particular, for Pocket 1, the average RMSD values of 2.27, 3.09, 2.58, and 2.59 Å ( Figure 8A) were computed for 14, 16, 83, and 85, respectively. On the other hand, for Pocket 2, we observed average RMSD values in the range of 2.17-2.51 Å ( Figure 8B). contacts shown as carbon sticks . Compounds 14, 16, 83, and 85 are depicted as pink, green, violet, and cyan carbon sticks, whereas H-bonds and π-π interactions are indicated as yellow and cyan dashed lines, respectively.
The best docked poses of 14, 16, 83, and 85 in complex with the most representative cluster of gD-Pocket 1 and gD-Pocket 2 were submitted to 500 ns of MDs using Desmond [44]. The results of MDs were investigated in terms of stability and conformational flexibility in the presence of the selected compounds. The stability of the complexes was evaluated by calculating the RMSD of the protein's backbone atoms from its initial to final conformation over the whole simulation. By RMSD analysis, we observed that the most promising compounds maintained overall stability throughout MDs in both Pockets, as shown in Figure 8. In particular, for Pocket 1, the average RMSD values of 2.27, 3.09, 2.58, and 2.59 Å ( Figure 8A) were computed for 14, 16, 83, and 85, respectively. On the other hand, for Pocket 2, we observed average RMSD values in the range of 2.17-2.51 Å ( Figure 8B). 63 kcal/mol, respectively. ADME parameters of the most promising compounds were predicted using the SwissADME server, and the obtained data are reported in Table 5. Furthermore, for each system, the binding free energy ∆Gbind was calculated using MM/GBSA methodology, extracting 100 snapshots from 500 ns of MDs. MM/GBSA analysis showed that the average calculated ∆Gbind of 14, 16, 83, and 85 complexed with gD-Pocket 1 were −38.88, −36.52, −40.22, and −31.14 kcal/mol, respectively, during the entire trajectories. The average ∆Gbind of 14, 16, 83, and 85 in complex to gD-Pocket 2 were −31.78, −29.71, −36.28, and −34.63 kcal/mol, respectively. ADME parameters of the most promising compounds were predicted using the SwissADME server, and the obtained data are reported in Table 5. All the investigated compounds are achiral, fit Lipinski's rule of five, and were not found to be pan assay interference compounds (PAINS), thus resulting suitable for further HSV-1 in vitro investigation.

MDs, GBPM, and MM/GBSA Calculations of gD-HVEM and gD-Nectin-1 Complexes
The best docked poses of each complex were submitted to 100 ns of MDs using Desmond ver. 4.2 [44]. To perform simulations in an aqueous biological environment, an appropriate system was built using OPLS_2005 as force field and an orthorhombic box with TIP4P water model extending of 10 Å outside the complex in all sides. The systems were maintained at a salt concentration of 0.15 M by adding appropriate Cl − counter ions to neutralize them to maintain the physiological condition. After optimization of the solvated models, we relaxed the systems with the Martyna−Tobias−Klein isobaric−isothermal ensemble (MTK_NPT). Finally, 100 ns unconstrained MDs were carried out using the following conditions: the NPT ensemble, a constant temperature of 300 K, a pressure of 1 bar, and a recording interval equal to 100 ps both for energy and for trajectory collecting 1000 frames for each simulation.
For both complexes, all frames were considered for GBPM analysis [38]. As previously reported [50], in order to evidence hydrophobic and hydrogen bond donors and acceptors spots, we used DRY, N1, and O GRID probes, respectively. For each complex, gD was seen as guest and HVEM and Nectin-1 as hosts. The selected residues at the interface covered a maximum distance of 3 Å from the most relevant interaction energy points (GBPM features) of the computed molecular interaction fields (MIFs). After selecting an energy cutoff 30% above the global minimum, the pivotal hotpots were resulted by the summa of the related GBPM features interaction energy.
One thousand snapshots from 100 ns of MDs were applied for the Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) free energy calculations [51] based on the following equation: where G comp , G pro and G lig denotes the free energy of the complex, protein, and the ligand; by splitting the energy contribution, it referred to ∆E ele , ∆E vdw and ∆E int as the gas-phase interaction energy between protein and ligand, thus including the electrostatic energy term, the van der Waals energy term, and the bond, angle, and dihedral terms, respectively. On the other hand, ∆E GB and ∆E surf indicate the polar and nonpolar desolvation free energy, respectively. The implicit solvation was calculated using the GB model [52], and the non-polar solvation energy was calculated using the solvent accessible surface area algorithm. The ∆Gbind reported in this study omitted the entropy contribution due to its relatively high computational demand and the lack of information of the conformational entropy that could lead to the introduction of additional error into the results [53]. The MDs trajectories were clustered based on the RMSD matrix of backbone atoms, and we obtained three representative structures for each complex. After removing the HVEM and Nectin-1 structures, we used a total of six gD conformations for the SBVS. The target binding sites were defined by a regular grid of about 20 Å centered on the residues responsible for binding with the cell receptor [8,40,41,54]. The residues that defined the gD binding site at the interface with HVEM, called for clarity "gD-Pocket 1", were as follows: Ala7, Ser8, Leu9, Lys10, Met11, Ala12, Asp13, Pro14, Asn15, Val24, Leu25, Asp26, Gln27, Leu28, Thr29, Asp30, Pro31, and Pro32. The gD binding site at the interface with Nectin-1, called "gD-Pocket 2", was characterized by the following residues: Pro23, Leu25, Gln27, Arg36, Val37, Tyr38, His39, Gln132, Val214, Asp215, Ser216, Ile217, Gly218, Met219, Leu220, Pro221, Arg222, Phe223, Thr230, Val231, and Tyr234.
The selected ligands were taken from an in-house small library of [1,2,3] Table 1). The 2D structures, reported in Table 1, were drawn using the ChemDraw Ultra 7.0 software and converted into 3D form using the import structures panel from Schrodinger maestro interface.
All compounds were optimized via the Ligprep module, [55] considering their ionization state at pH 7.4, and energy minimized using OPLS_2005 as force field [37].
As reported in other studies [56][57][58][59], the docking simulations of our focused library were computed using the Glide [60] ligand flexible algorithm at the standard-precision (SP) level, generating 10 possible poses for each site. The best docked poses for gD-Pocket 1 and gD-Pocket 2 were submitted to 500 ns of MDs in order to define the structural and energy profile of the best ligands in complex with both gD-Pockets. The simulations were carried out under the above-mentioned conditions. All simulations were performed by Desmond package [44] and "Simulation Interactions Diagram" panel was used as a post-MD analysis tool for exploring protein-ligand interactions. MM/GBSA free energy calculations of the best generated complexes were conducted along 100 frames of 500 ns of MDs. ADME descriptors and pharmacokinetic properties of the promising compounds were predicted by means of SwissADME tool [61].

Conclusions
In the present computational investigation, the pivotal interactions occurring at the interface of gD and its human surface receptors were carefully explored. Since the avail-ability of X-ray structures of gD in complex with HVEM and Nectin-1, this study aimed at highlighting the key residues involved in the binding interface during MDs. Our reasoning was based on the idea that molecular recognition entails a two-way influence between the interacting partners, whereby the flexibility of gD and its human receptors determined the optimal conformation for the complexes formation. This aspect was better appreciated analyzing the behavior of the complexes after MDs, and the resulting dynamic adaptation was used to define the principal residues responsible for forming a stable complex. MM/GBSA analysis revealed a greater binding affinity of gD towards Nectin-1 with respect to HVEM, given that the rearrangement of N-terminal hairpin is not necessary for gD-Nectin-1 interaction, which instead ensured a more rigid gD-binding pocket. By analyzing HSV-1 PPIs, through the GBPM method, gD N-terminus residues at position 10,12,13,14,15,26,27,29,33, and 35 were found pivotal for gD-HVEM binding, whereas nine residues, at position 26, 27, 38, 215, 218, 219, 220, 221, and 222 were relevant for the connection with Nectin-1. Afterwards, an SBVS targeting different conformations of HSV-1 gD binding interfaces was applied. Developing PPI inhibitors is challenging [62], owing to issues such as the general lack of small-molecule starting points for drug design, the typical flatness of the interface, the difficulty of distinguishing real from artefact binding, and the size and character of traditional libraries. In silico approaches are a consolidated strategy to speed up the drug discovery process, as demonstrated during the pandemic emergency, increasing our understanding of how biological systems work and translating this knowledge into new molecules with interesting therapeutic potential. In this context, the inhibition of HSV-1 through targeting viral gD protein represents a good way to act with virus adsorption and membrane fusion. Herein we identified four triazole-pyridines as potential HSV-1 gD inhibitors, with a favourable pharmacokinetic profile and a good theoretical affinity towards all conformations of HSV-1 gD. The surface plasmon resonance technique could provide further information into the mechanism of action at a molecular level and could move forward novel antiviral compounds. Our study suggests a promising basis for the design of a new generation of anti-HSV-1 drugs targeting gD-receptor interfaces.