Edit · 2026 01 10 132856 Nature S41598 025 87407 8

Saving needs a token. Log in or register to get one. Edits are public and licensed CC-BY-SA-4.0.

# Unbiased picture of the ligand docking process for the hevein protein–oligosaccharide complex | Scientific Reports

**URL:** https://www.nature.com/articles/s41598-025-87407-8

**Saved:** 2026-01-10T21:29:01.868Z

---

Download PDF
Article
Open access
Published: 27 January 2025
Unbiased picture of the ligand docking process for the hevein protein–oligosaccharide complex
Toshifumi Yui & Takuya Uto

Scientific Reports 
volume
 15, Article number: 3335 (2025) Cite this article

2101 Accesses

Metrics
details

Abstract

The ligand-docking behavior of hevein, the major latex protein from the rubber tree Hevea brasiliensis (Euphorbiaceae), has been investigated by the unguided molecular dynamics (MD) simulation method. An oligosaccharide molecule, initially placed in an arbitrary position, was allowed to move around hevein for a prolonged simulation time, on the order of microseconds, with the expectation of spontaneous ligand docking of the oligosaccharide molecule to the binding site of hevein. In the binary solution system consisting of a hevein molecule and a chito-trisaccharide (GlcNAc3) molecule, three out of the six separate simulation runs successfully reproduced the complex structure of the observed binding from. It appeared that the surface topology formed by two aromatic side chains of the hevein molecule played a role in orienting the GlcNAc3 molecule in the correct direction. We also performed MD simulations of the ternary solution system containing a cello-hexasaccharide (Glc6) molecule in addition to hevein and a chito-hexasaccharide (GlcNAc6) molecule. Formation of hevein–GlcNAc6 complex structures was exclusively observed, while the Glc6 molecule remained in the solvent phase throughout the simulations. Obviously, the acetamide groups of GlcNAc play a role in detecting the binding site and its vicinity on the protein surface.

Similar content being viewed by others
Comprehensive structural assignment of glycosaminoglycan oligo- and polysaccharides by protein nanopore
Article Open access
30 August 2022
Molecular mechanism of decision-making in glycosaminoglycan biosynthesis
Article Open access
13 October 2023
Targeting protein–ligand neosurfaces with a generalizable deep learning tool
Article Open access
15 January 2025
Introduction

Many plants protect themselves from pathogenic attack by producing defense proteins, some of which have the ability to reversibly bind to chitin, β-(1→4)-linked N-acetyl D-glucosamine (GlcNAc), which functions as a structural polysaccharide in the cell walls of fungi and the exoskeletons of invertebrates. Most plant proteins with antifungal activity contain a noncatalytic, plant-specific chitin-binding domain, which is called the hevein domain after hevein, the major latex protein of the rubber tree, Hevea Brasiliensis (Euphorbiaceae), involved in the coagglutination of rubber particles1. This domain shares a common structural motif of 30–43 residues rich in glycine and cysteine residues at highly conserved positions and organized around a three to five conserved disulfide bond core2. In addition to hevein itself, this domain has been detected in pseudohevein, Urtica dioica agglutin, wheat germ agglutin, Amaranthus caudatus antimicrobial peptides, and the chitin-binding modules of class I chitinases3,4.

Three-dimensional structure analyses of hevein and its chito-oligosaccharide complexes have been performed by NMR and X-ray measurements, which have provided the atomistic details of the ligand-binding structures5,6,7,8,9. Aboitiz et al.10 proposed the three-dimensional NMR structure of HEV32, truncated hevein lacking eleven C-terminal amino acids, and the triacetylchitotriose (GlcNAc3) complex form. In the structure of the hevein–GlcNAc3 complex (Fig. 1), two of the tryptophan side chains of the hevein molecule (W21 and W23) are bound to the pyranose rings of the GlcNAc3 molecule, which appears to be stabilized by CH–π stacking interactions and van der Waals interactions at the interface of the complex. The GlcNAc residue at the nonreducing end is attached to the side chain of W23, which can be defined as the subsite + 1. The hydroxyl groups of residues serine 19 (S19) and tyrosine 30 (Y30) are involved in hydrogen bonding with the carbonyl group of the acetamide moiety and the O3–H group of the nonreducing GlcNAc residue.

Fig. 1

Reference model of the complex of hevein with a chito-trisaccharide (GlcNAc3) molecule that is accommodated at subsites + 1 to + 3 from a nonreducing residue10. The three aromatic amino acids composing the binding site are shown in space-filling representation.

Full size image

Many noncatalytic carbohydrate-binding proteins, including hevein, provide a shallow ligand-binding surface11 for carbohydrate ligands whose polar functional groups are almost symmetrically aligned with respect to the molecular chain, suggesting that the proteins can form different ligand binding modes with similar binding energies. Jiménez-Barbero and co-workers12,13 discussed the docking characteristics of the hevein complex, which they referred to as a typical example of protein–carbohydrate interactions. In particular, they pointed out that this feature can be attributed to the amphiphilic character of carbohydrate ligands12,13. The hydroxyl groups with hydrophilic nature and the aliphatic C–H region forming a hydrophobic patch on the pyranose ring faces account for the formation of hydrogen bonds with the side-chain of the polar amino acid and the occurrence of hydrophobic interactions with the aromatic side chain, respectively. The latter interactions result from the mutual shielding from bulk water at the interface. The interfaces may also involve attractive electrostatic interactions derived from the positive net charge of the C–H groups and the delocalized π electrons of the aromatic rings.

Theoretical calculations, such as molecular dynamics (MD) simulations and quantum mechanics calculations, have provided further insights into the motional behavior at the binding site, including the bound oligosaccharide ligand conformation and thermodynamics associated with ligand binding10,12,14,15,16,17,18. These theoretical studies were performed based on the experimental structure as a starting structure or as a target structure to simulate the process of ligand docking. Shan et al.19 reported an unguided MD simulation study of the prediction of the complex structure of Src kinase and the kinase inhibitor PPI, in which the MD simulations were continued until the PPI ligand initially located at an arbitrary position in the simulation box reached the target binding site on the protein surface. Several of the simulations successfully provided the ligand-binding structure, which was virtually identical to the crystallographically determined complex structure. The simulation trajectories provided a continuous atomistic view of the entire ligand-docking process, including the intermediate ligand conformations and the role of solvent water molecules in the ligand binding. The same simulation tactic was applied to the MD system consisting of HEV32 and chito-oligosaccharide molecules, either GlcNAc, GlcNAc2, or GclNAc3, by Solanke et al.20 While most of the simulation trajectories showed irreversible binding of the chito-oligosaccharide molecules to the ligand binding sites, the MD simulation with a GlcNAc3 molecule failed to reproduce the identical ligand-binding mode to that observed in the NMR complex structure10. Instead, the reducing-end residue of the GlcNAc3 molecule was bound to the subsite + 1 while the disaccharide moiety remained in the bulk water.

In the present study, we performed unguided MD simulations combined with modified van der Waals parameters and extended the simulation time to 10 µs. We first confirmed that the MD simulations successfully reproduced the identical hevein–GlcNAc3 complex structure to that observed in the NMR complex structure10. The ligand-docking processes were analyzed by visualizing the MD trajectories. In particular, we investigated the possibility that a GlcNAc3 molecule binds to the substrate-binding site via a specific pathway during complex formation, including the possibility of the presence of a hot spot; this information cannot be obtained by means of experimental approaches such as NMR and X-ray measurements. The MD simulations were then applied to predict the structures of an unknown complex form with a longer substrate, a hevein–GlcNAc6 complex, and to evaluate the substrate specificity of hevein to a chito-oligosaccharide compared to a cello-oligosaccharide.

Methods

The X-ray crystal structure of a hevein domain (4WP4) extracted from latex surgical gloves9 was used as the starting structure. The protein consisted of 43 amino acids and involved four disulfide bonds. The reference model of the carbohydrate–hevein complex was constructed by combining the two structures of a Hev b 6.02 molecule and a truncated hevein (HEV32) complex with a chito-trisaccharide (GlcNAc3) (1T0W)10. The binding site of hevein is composed of the aromatic side chains W23 and W21, which serve as the subsite + 1 and + 2 for a chitin substrate, respectively (Fig. 1). The + 1 subsite accommodates the nonreducing GlcNAc residue, and its 2-acetamide group is the contact with the side chain of Y30. Two hydrogen bonds are formed between the 3-OH group of the nonreducing residue and the hydroxyl groups of the Y30 and S19 side chains10. Six binary solution systems consisting of a hevein molecule and an oligosaccharide ligand were constructed by arbitrarily placing an oligosaccharide ligand on each of the six sides of the hevein molecule (Fig. 2). The oligosaccharides used in the present MD simulations were a chito-trisaccharide (GlcNAc3), a chito-hexasaccharide (GlcNAc6), and a cello-hexasaccharide (Glc6). The latter two saccharides were used in competitive unguided MD simulations for ternary solution systems in which GlcNAc6 and Glc6 molecules were first placed on opposite sides of a hevein molecule. These solute systems were solvated with 5000 TIP3P water models21,22.

The solvated complex systems were subjected to energy minimization steps in which the water configurations were first optimized, followed by full structure optimizations. The systems were gradually heated from 100 to 300 K in 500 ps using the NVT ensemble. The density of the systems was equilibrated by 500-ps NPT ensemble simulations with a constant temperature of 300 K and pressure of 1 bar. The equilibrated systems were subjected to unbiased production MD simulations using the NPT ensemble for 10 and 5 µs for the hevein–GlcNAc3 solution, and hevein–GlcNAc6 and hevein–GlcNAc6/Glc6 solutions, respectively. For both the binary and ternary solution systems, six separate simulations were performed, each differing in the initial position(s) of the oligosaccharide ligand(s). The temperature and pressure of the system were regulated by a Langevin thermostat23,24 with a collision frequency of 5 ps− 1 and a Berendsen barostat25 with a pressure relaxation time of 1.0 ps, respectively. A time step of 2.0 fs was used with the SHAKE option26 to ensure rigid hydrogen atoms. The particle mesh Ewald method27 was used for the long-range interactions, and the van der Waals interactions were cut off at 8.0 Å. All of the energy minimizations and MD simulations were performed with the PMEMD and PMEMD.CUDA modules28,29 of the Amber 16 program package30 using an NVIDIA Kepler GPU system. The hevein molecule was described by the ff14SB force field31,32,33, and the oligosaccharide molecules were described by the GLYCAM06 force field34. In the present MD simulations, the parameters describing the intermolecular interactions of the carbohydrate and amino group carbon atoms were corrected to weaken the excessively attractive interactions by using modified Lennard-Jones parameter values developed to reproduce the experimental osmotic pressures of mixed aqueous solutions of diglycine and sucrose35,36.

The spatial motion of the oligosaccharide ligand(s) during the MD simulations was monitored by evaluating the root-mean-square deviation (RMSd) values of the positions of the pyranose ring and glycosidic oxygen atoms, with respect to those of the observed complex NMR structure10 as a reference structure. The value indicated how far an oligosaccharide ligand was located from the binding site, so that when the bound position of a GlcNAc3 part was identical to that in the NMR structure RMSd = 0 Å. As an alternative assessment of the bound position, the binding energy (Ebind), representing the nonbonded interactions composed of van der Waals and electrostatic interactions between an oligosaccharide molecule and the hevein surface in the gas phase, and the binding free energy including Poisson-Boltzmann solvation energy of an implicit solvent model (GbindPBSA) were calculated from the MD trajectories using the MMPBSA module37. The MD trajectories were analyzed using the CPPTRAJ module of AmberTools 1938, and VMD 1.9.339 was used for visual analysis. The molecular graphics software PyMOL (version 2.5.0)40 was used to visualize the system.

Fig. 2

Six starting positions of the oligosaccharide (GlcNAc3) molecule around the hevein molecule.

Full size image
Results and discussion
Ligand-docking behavior of the binary solution containing hevein and GlcNAc3 molecules

The ligand-docking behavior in the MD simulations of the binary solution containing hevein and GlcNAc3 molecules is shown in Fig. 3, which is represented by the two-dimensional ligand-RMSd and binding energy (ΔEbind) distributions (Fig. 3A) and ligand-RMSd profiles with respect to the simulation time (Fig. 3B). Among the six separate docking simulations, simulation runs 2, 3, and 6 showed the formation of the hevein–GlcNAc3 complex, while hevein remained in the unbound state in the remaining simulation runs 1, 4, and 5 throughout 10 µs simulation times. In Fig. 3A, the three high-density dots with RMSd values close to the origin, namely, those near the reference complex structure, indicate the formation of the complex structures. All of the complex structures, labeled bound poses a, b, and c, survived for approximately 100–400 ns before they reversibly changed to an unbound state (Fig. 3B). The two complex structures of bound poses a and b obtained from the MD trajectories of simulation runs 2 and 3 are shown in Fig. 3C. Bound pose a, in which the GlcNAc residues were accommodated at subsites + 1 to + 3, was consistent with that of the observed NMR complex structure10. Although the complex structures of bound pose a appeared in higher ΔEbind ranges than those of bound poses b and c, bound pose a can be considered to be the most stable in terms of the total complex lifetime. It should be noted that the contribution of the solvent effect was not taken into account in the calculation of ΔEbind, which is the sum of the intermolecular van der Waals and electrostatic interactions in the gas phase. Compared with bound pose a, the nonreducing end residue of the GlcNAc3 molecule shifted to subsite − 1 from subsite + 1 in bound pose b. A small distribution of ligand-RMSd and ΔEbind dots, labeled c, appeared at approximately 8500 ns for simulation run 2 with a relatively short lifetime. Bound pose c was identical to bound pose b with respect to the subsites occupied by the GlcNAc residues, but the direction of the GlcNAc3 molecule was reversed.

As mentioned previously, the structure of hevein is distinguished by rich in conserved cysteine and glycine residues, which are a common structural motif observed in hevein domains2. Hevein involves five disulfide bonds to secure stability of the three-dimensional structure, whereas glycine residues may contribute flexibility to the peptide backbone. Fig. S1 shows a hevein backbone-RMSd profile with respect to the simulation time obtained from the final 5000 ns MD trajectory of the simulation run 3. The observed complex NMR structure10 was used as a reference structure. To emphasize RMSd variations, the residues subject to the RMSd calculations were those belonging to the amino acid sequence between S19 and Y30. When comparing the RMSd profiles between the state of the bound pose a and an unbound state, a slight reduction in RMSd values and the absence of distinct pulse peaks are observed in the bound state. About 200 ns after the time of the complex formation of bound pose b, the RMSd values increased discontinuously by about 0.25 Å. This indicates that the excessive binding energy occurring in the complex structure of bound pose b led to local steric distortions of the peptide backbone, which could increase the whole steric energy of the hevein structure.

The three ligand docking processes leading to the complex structures of bound pose a were visually examined by inspecting the MD trajectories taken from simulation runs 2, 3, and 6. In the trajectory from 657 to 664 ns of simulation run 2 (Movie S1), the GlcNAc3 molecule showed a multistage process as it bound to the binding site. When it first approached the binding site, the GlcNAc3 molecule appeared to climb up the side chain of Y30 as if it were a handhold. The GlcNAc3 molecule first attached to the side chain of W21 and then made a combined movement of flipping and turning on the aromatic side chains to fit into the binding site with bound pose a. In the trajectory from 3316 to 3325 ns of simulation run 6, the GlcNAc3 molecule directly reached the binding site and remained on the upper edge of the aromatic side chains, where the aromatic rings slightly overlapped for a few nanoseconds. The GlcNAc3 molecule then docked into the binding site in concert with the elimination of the overlapping part to better interact with the binding surface (Movie S2). The ligand docking process in the trajectory from 9677 to 9695 ns of simulation run 3 was similar to that observed in simulation run 6. The GlcNAc3 molecule bound to the binding site from the top edge of the aromatic side chains in a rather straightforward manner (Movie S3). The ligand docking motion shown in the three movies suggests that the surface topology formed by the two aromatic side chains W21 and W23 primarily regulates a GlcNAc3 molecule to bind in a proper fashion, especially in terms of its direction with respect to the order of the subsites. The ligand docking motion shown in Movies M2 and M3 suggested that the intermolecular hydrogen bonds involving the hydroxyl groups of the S19 and Y30 side chains and the acetamide moiety of the nonreducing end residue may play a role in reinforcing the complex structure after the GlcNAc3 molecule correctly binds to the binding site.

Our initial unguided MD simulation study of the hevein–GlcNAc3 binary solution system using the original nonbonded interaction parameters was affected by the excessively high affinity between the sugar and amino acid residues. In all of the six simulation runs, the GlcNAc3 molecule irreversibly bound to the binding site in the course of 1 µs simulations only to provide the undesired bound pose b with RMSd values of approximately 5 Å (see Fig. S2). From Fig. 3C, the ligand binding structure of bound pose b appears to involve a larger binding interface than that of bound pose a, which is consistent with lower ΔEbind values of the complex structures of bound pose b, regardless of the nonbonded interaction parameters used. Clearly, the solvent effect is an additional effect that makes the complex structure of bound pose a the most stable. Table S1 lists the thermodynamic values calculated from the MD trajectory frames involving the bound poses a and b complex structures. The average values of ΔEbind for bound pose a and b calculated from the corresponding trajectory frames of simulation runs 2, 3, and 6 were − 11.8 ± 2.4 kcal/mol and − 26.5 ± 4.0 kcal/mol, respectively. Aboitiz et al. reported the van’t Hoff enthalpy of -8.7 kcal/mol for the hevein–GlcNAc3 complex formation determined by the NMR titration measurement, in addition to the observed binding free energy of -5.5 kcal/mol10. However, the ΔGbindPBSA calculations provided undesirable positive values, suggesting that the modified van der Waals parameters may not have worked properly on the MM_PBSA.py module. On the other hand, the MD simulations using the original van der Waals parameters for the NMR complex structure suggested an excessively low ΔEbind value of -51.4 kcal/mol, accompanied by the apparently reasonable ΔGbindPBSA value of -22.8 kcal/mol. It should be noted that in the present MD simulations the evaluation of the thermodynamic values, in combination with the RMSd evaluations of the ligand positions, was rather adopted as a means of monitoring complex formation.

Fig. 3

Ligand docking behavior of the binary solution containing hevein and GlcNAc3 molecules obtained from the trajectories of simulation runs 2, 3, and 6. (A) Two-dimensional ligand-RMSd (Å) and binding energy (ΔEbind, kcal/mol) distributions. Each dot represents a spatial arrangement of the two solute molecules with the average ligand-RMSd and ΔEbind values obtained from a 1 ns MD trajectory. The labels a, b, and c denote the bound poses. (B) Ligand-RMSd profiles with respect to the simulation time (ns). (C) Representative complex structures of bound poses a and b.

Full size image
Ligand-docking behavior of the binary solution containing hevein and GlcNAc6 molecules

Originally, hevein binds to the macromolecular chitin that composes the fungal cell wall. Therefore, the question arises as to whether the ligand docking behavior observed for the GlcNAc3 binary solutions is inherently different from that for polymeric substrates. Unguided MD simulations of a binary solution containing a hevein molecule and a GlcNAc6 molecule were performed in a similar fashion to obtain a more realistic picture of the ligand docking behavior. Compared with the docking behavior observed for the hevein–GlcNAc3 binary solutions, the GlcNAc6 molecule more effectively bound to the binding site in four out of the six simulation runs during a simulation time of 5 µs, as indicated by the two-dimensional ligand-RMSd and ΔEbind distributions (Fig. 4A) and the time course of the ligand-RMSd value (Fig. 4B). Three bound poses were detected, all showing an identical direction of the GlcNAc6 molecule to that of the observed hevein–GlcNAc3 complex structure (Fig. 4C). All of the three subsites, + 1 to + 3, were occupied by GlcNAc residues in bound pose a. Note that the RMSd values, using the observed hevein–GlcNAc3 complex as the reference structure, reflect the deviations in the atomic coordinates of the GlcNAc residues accommodated at the three subsites from those in the observed structure. The bound poses with one vacant subsite, bound poses b and c, tended to have relatively large RMSd values. Similar to the ligand-docking behavior observed for the hevein–GlcNAc3 binary solutions, neither transient ligand binding to the protein surface other than to the binding site nor a specific intermediate complex structure was observed. It appeared that the complex structures of bound pose b were the most stable, as suggested by their long lifetimes of more than 99% and 60% of the simulation time in simulation runs 2 and 5, respectively, accompanied by ΔEbind values of mostly approximately − 40 kcal/mol. The GlcNAc6 molecule unevenly bound to the binding site in all three bound poses with either the reducing or nonreducing moiety of the GlcNAc6 molecule exposed to the bulk water (Fig. 4C). It is likely that the solvent effect acting on the unbound three to four residue moiety of the GlcNAc6 molecule promotes complex formation.

Fig. 4

Ligand docking behavior of the binary solution containing hevein and GlcNAc6 molecules obtained from the MD trajectories of simulation runs 2, 3, 4, and 5. (A) Two-dimensional ligand-RMSd (Å) and binding energy (ΔEbind, kcal/mol) distributions. Each dot represents a spatial arrangement of the two solute molecules with the average ligand-RMSd and ΔEbind values obtained from a 1 ns MD trajectory. The labels a, b, and c denote the bound poses. (B) Ligand-RMSd profiles with respect to the simulation time (ns). (C) Representative structures of bound poses a, b, and c.

Full size image
Competitive ligand-docking behavior of a ternary solution containing hevein, GlcNAc6, and Glc6 molecules

The present unguided MD simulation method for a binary solution system containing a ligand molecule and its target protein was originally proposed for identifying an unknown binding site and predicting a complex structure19. As an extended application of the unguided MD simulation method, we propose using a ternary solution system containing hevein, an oligosaccharide substrate, and a substrate analog to compare their relative affinities for the binding site, where the two ligand molecules are expected to competitively bind to the binding site. The oligosaccharide ligands used were GlcNAc6 and Glc6 molecules, the latter playing the role of a substrate analog molecule. A hevein–GlcNAc6 complex successfully formed in two out of the six simulation runs, while the Glc6 molecule remained circling around the hevein molecule throughout the 5 µs simulation time in all of the six simulation runs (Fig. 5A and B). In the two hevein–GlcNAc6 complex structures, all three subsites were occupied by GlcNAc residues (Fig. 5C). The bound pose b corresponds to the bound pose a obtained in the previous MD simulations of the GlcNAc6 binary solution system. The position of the GlcNAc residue at the subsite + 3 slightly deviated from that of the corresponding residue in the reference complex structure, resulting in larger RMSd values. The nonreducing residue of the GlcNAc6 molecule shifted from the subsite + 1 to the subsite − 1 by changing the ligand-binding mode from bound pose b to bound pose a. We did not anticipate that the competitive MD simulations would be able to distinguish the two oligosaccharide molecules so clearly, owing to the limited simulation system. Interestingly, the complex structures with the GlcNAc6 molecules were retained for longer simulation times (up to a few microseconds) than those in the GlcNAc6 binary solutions (hundreds of nanoseconds), which was obviously caused by the double effective concentration of the GlcNAc6 molecule in the ternary solution. It should be noted that both of the oligosaccharide molecules with the common β-1→4 glycosidic linkage prefer to have a similar extended overall structure, even in the solution state. Therefore, the competitive MD simulations again confirmed that the 2-acetamide groups significantly contribute to the substrate specificity of the hevein binding site. This group may act as a guide to lead a chito-oligosaccharide to the hevein binding site, in addition to forming intermolecular hydrogen bonds with the hydroxyl groups of the S19 and Y30 side chains10, probably mediated by water molecules.

Fig. 5

Competitive docking behavior of the ternary solution containing hevein, GlcNAc6, and Glc6 molecules obtained from the MD trajectories of simulation runs 2 and 3. (A) Two-dimensional ligand-RMSd (Å) and binding energy (ΔEbind, kcal/mol) distributions. Each dot represents a spatial arrangement of the two solute molecules, hevein and GlcNAc6, with the average ligand-RMSd and ΔEbind values obtained from a 1 ns MD trajectory. The labels a and b denote the bound poses of the hevein and GlcNAc6 complex structures. (B) Ligand-RMSd profiles with respect to the simulation time (ns). (C) Representative structures of bound poses a and b, accompanied by a Glc6 molecule.

Full size image
Conclusions

The present unguided MD simulations of a binary solution system containing GlcNAc3 and hevein molecules successfully predicted the complex structure of the observed bound pose10 as the most stable form. The MD simulations also gave alternative bound pose b in which the GlcNAc3 molecule, shifted by one residue, was accommodated in the subsites − 1 to + 2 with the same orientation. Using weakened nonbonded interactions with modified van der Waals parameters, combined with the solvent effect, allowed the GlcNAc3 molecule to reversibly bind to the binding site for several hundred nanoseconds. In the hevein–GlcNAc6 binary solution systems, the GlcNAc6 molecule more effectively but unevenly bound to the binding site. The terminal moiety of the GlcNAc6 molecule was exposed to the bulk water, suggesting that the nonbonded interactions at the binding interface and the solvent effect on the unbound GlcNAc residues were balanced to form a stable complex structure. In our previous docking MD simulations using the original van der Waals parameters, in some simulation runs, a chito-oligosaccharide molecule bound to the hevein surface other than to the binding site. Such an undesirable complex structure was not observed in the present MD simulations, indicating that using the modified van der Waals parameters allowed the chito-oligosaccharide molecule to successfully identify the binding site from the rest of the protein surface.

The docking processes of a chito-oligosaccharide molecule were evaluated by visualizing the MD trajectories. After moving around the hevein molecule, the chito-oligosaccharide molecule rapidly approached the binding site from a different direction in each simulation run, without a specific hot spot where the ligand molecule was temporarily bound before moving to the binding site. It appeared that the chito-oligosaccharide molecule approaching the binding site was oriented in the appropriate bound pose on its own accord, although it was in fact mainly guided by the surface topology of the aromatic side chains W21 and W23. While intermolecular hydrogen bonds with the side chains of S19 and Y30 are then formed to secure the complex structure, as mentioned above, surrounding water may assist in guiding the GlcNAc3 molecule to the binding site and stabilize the final complex structure.

We also performed competitive docking simulations for a ternary solution containing hevein and two different oligosaccharide substrates, GlcNAc6 and Glc6 molecules, as an extended application of the unguided MD simulation method. To our surprise, hevein clearly distinguished the GlcNAc6 molecule from the Glc6 molecule, yielding exclusively a hevein–GlcNAc6 complex. In addition to participating in the formation of the intermolecular hydrogen bonds, the acetamide groups, which are bulkier than the hydroxyl groups of a glucose residue, may act as a detector to explore the binding site and its vicinity on the protein surface. However, we are reluctant to conclude that the binding site of hevein has no affinity for the Glc6 molecule at all, considering that the present MD simulations were performed for limited simulation time on the order of nanoseconds or microseconds, a short time compared with the real time scale. One of the drawbacks of the unguided MD simulation method is that it is essentially impossible to rule out complex formation for a given ligand. When performing a similar competitive docking simulation for hybrid substrates with closer affinities for a target protein, the simulation conditions, such as the solute concentration and simulation time, must be carefully examined; otherwise, a false conclusion could be drawn from an accidental result.

The present docking MD simulations were not intended to be an unbiased search for a ligand binding site of a target protein, contrary to the original objective of the unguided MD simulation method19. Because hevein is a small protein and its binding site has been defined, the purpose of the present MD simulation study was to confirm the reproducibility of the observed hevein–GlcNAc3 complex structure, which was successfully achieved by using modified van der Waals parameters to weaken the nonbonding interactions between the GlcNAc3 and hevein molecules. With this information, it is expected that a future MD study of the hevein complex using an advanced MD method, such as thermodynamic integration, will be performed with high reliability.

Data availability

The video files, Movie S1, S2, and S3, and the word file containing Fig. S1 and S2, and Table S1 are available in the Supplementary Material. The AMBER topology, input, and initial coordinates, and the PDB files of the starting structures used for the present study are also provided in the Supplementary Material.

References

Gidrol, X., Chrestin, H., Tan, H. L. & Kush, A. Hevein, a lectin-like protein from Hevea brasiliensis (rubber tree) is involved in the coagulation of latex. J. Biol. Chem. 269, 9278–9283 (1994).

Article
 
CAS
 
PubMed
 
Google Scholar

Drenth, J., Low, B. W., Richardson, J. S. & Wright, C. S. The toxin-agglutinin Fold. A new group of small protein structures organized around a four-disulfide core. J. Biol. Chem. 255, 2652–2655 (1980).

Article
 
CAS
 
PubMed
 
Google Scholar

Damme, E. J. M. V., Peumans, W. J., Barre, A. & Rougé, P. Plant lectins: a composite of several distinct families of structurally and evolutionary related proteins with diverse biological roles. Crit. Rev. Plant. Sci. 17, 575–692 (1998).

Article
 
Google Scholar

Beintema, J. J. Structural features of plant chitinases and chitin-binding proteins. FEBS Lett. 350, 159–163 (1994).

Article
 
CAS
 
PubMed
 
MATH
 
Google Scholar

Andersen, N. H., Cao, B., Rodriguez-Romero, A., Arreguin, B. & Hevein NMR assignment and assessment of solution-state folding for the agglutinin-toxin motif. Biochemistry 32, 1407–1422 (1993).

Article
 
CAS
 
PubMed
 
Google Scholar

Asensio, J. L., Canada, F. J., Bruix, M., Rodriguez-Romero, A. & Jimenez-Barbero, J. The Interaction of Hevein with N-acetylglucosamine-containing oligosaccharides. Solution structure of hevein complexed to chitobiose. Eur. J. Biochem. 230, 621–633 (1995).

CAS
 
PubMed
 
Google Scholar

Rodriguez-Romero, A., Ravichandran, K. G. & Soriano-Garcia, M. Crystal structure of hevein at 2.8 Å resolution. FEBS Lett. 291, 307–309 (1991).

Article
 
CAS
 
PubMed
 
MATH
 
Google Scholar

Reyes-Lopez, C. A. et al. Insights into a conformational epitope of hev b 6.02 (hevein). Biochem. Biophys. Res. Commun. 314, 123–130 (2004).

Article
 
CAS
 
PubMed
 
MATH
 
Google Scholar

Galicia, C., Mendoza-Hernandez, G. & Rodriguez-Romero, A. Impact of the vulcanization process on the structural characteristics and IgE recognition of two allergens, Hev b 2 and Hev b 6.02, extracted from latex surgical gloves. Mol. Immunol. 65, 250–258 (2015).

Article
 
CAS
 
PubMed
 
Google Scholar

Aboitiz, N. et al. NMR and modeling studies of protein-carbohydrate interactions: synthesis, three-dimensional structure, and recognition properties of a minimum hevein domain with binding affinity for chitooligosaccharides. ChemBioChem 5, 1245–1255 (2004).

Article
 
CAS
 
PubMed
 
MATH
 
Google Scholar

Boraston, A. B., Bolam, D. N., Gilbert, H. J. & Davies, G. J. Carbohydrate-binding modules: fine-tuning polysaccharide recognition. Biochem. J. 382, 769–981 (2004).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Jiménez-Barbero, J. et al. Protein-carbohydrate interactions: A combined theoretical and NMR experimental approach on carbohydrate-aromatic interactions and on pyranose ring distortion, in NMR Spectroscopy and Computer Modeling of Carbohydrates, ACS Symp. Ser. 930 (eds J. F. G Vliegenthar & R. J. Woods), 60–80 (2006).

Jiménez-Barbero, J. et al. Hevein domains: An attractive model to study carbohydrate–protein interactions at atomic resolution. In Adv. Carbohydr. Chem. Biochem., 60 (ed D. Horton), 303–354 (2006).

Colombo, G., Meli, M., Canada, J., Asensio, J. L. & Jimenez-Barbero, J. Toward the understanding of the structure and dynamics of protein-carbohydrate interactions: molecular dynamics studies of the complexes between hevein and oligosaccharidic ligands. Carbohydr. Res. 339, 985–994 (2004).

Article
 
CAS
 
PubMed
 
Google Scholar

Colombo, G., Meli, M., Canada, J., Asensio, J. L. & Jimenez-Barbero, J. A dynamic perspective on the molecular recognition of chitooligosaccharide ligands by hevein domains. Carbohydr. Res. 340, 1039–1049 (2005).

Article
 
CAS
 
PubMed
 
Google Scholar

Mareska, V., Tvaroska, I., Kralova, B. & Spiwok, V. Molecular simulations of hevein/(GlcNAc)3 complex with weakened OH/O and CH/π hydrogen bonds: implications for their role in complex stabilization. Carbohydr. Res. 408, 1–7 (2015).

Article
 
CAS
 
PubMed
 
Google Scholar

Porto, W. F., Souza, V. A. & Nolasco, D. O. Franco, O. L. In silico identification of novel hevein-like peptide precursors. Peptides 38, 127–136 (2012).

Article
 
CAS
 
PubMed
 
Google Scholar

Koppisetty, C. A., Frank, M., Lyubartsev, A. P. & Nyholm, P. G. Binding energy calculations for hevein-carbohydrate interactions using expanded ensemble molecular dynamics simulations. J. Comput. Aided Mol. Des. 29, 13–21 (2015).

Article
 
ADS
 
CAS
 
PubMed
 
Google Scholar

Shan, Y. et al. How does a drug molecule find its target binding site? J. Am. Chem. Soc. 133, 9181–9183 (2011).

Article
 
CAS
 
PubMed
 
PubMed Central
 
MATH
 
Google Scholar

Solanke, C. O. et al. Atomistic simulation of carbohydrate-protein complex formation: Hevein-32 domain. Sci. Rep. 9, 18918 (2019).

Article
 
ADS
 
CAS
 
PubMed
 
PubMed Central
 
MATH
 
Google Scholar

Jorgensen, W. L. Quantum and statistical mechanical studies of liquids. 10. Transferable intermolecular potential functions for water, alcohols, and ethers. Application to liquid water. J. Am. Chem. Soc. 103, 335–340 (1981).

Article
 
CAS
 
Google Scholar

Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W. & Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935 (1983).

Article
 
ADS
 
CAS
 
Google Scholar

Hoover, W. G., Ladd, A. J. C. & Moran, B. High-strain-rate plastic flow studied via nonequilibrium molecular dynamics. Phys. Rev. Lett. 48, 1818–1820 (1982).

Article
 
ADS
 
CAS
 
Google Scholar

Evans, D. J. Computer experiment for nonlinear thermodynamics of couette flow. J. Chem. Phys. 78, 3297–3302 (1983).

Article
 
ADS
 
CAS
 
MATH
 
Google Scholar

Berendsen, H. J. C., Postma, J. P. M., van Gunsteren, W. F., DiNola, A. & Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81, 3684–3690 (1984).

Article
 
ADS
 
CAS
 
Google Scholar

Ryckaert, J. P., Ciccotti, G. & Berendsen, H. J. C. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23, 327–341 (1977).

Article
 
ADS
 
CAS
 
MATH
 
Google Scholar

Essmann, U. et al. A smooth particle mesh Ewald method. J. Chem. Phys. 103, 8577–8593 (1995).

Article
 
ADS
 
CAS
 
MATH
 
Google Scholar

Le Grand, S., Götz, A. W. & Walker, R. C. Speed without compromise—A mixed precision model for GPU accelerated molecular dynamics simulations. Comput. Phys. Commun. 184 SPFP, 374–380 (2013).

Article
 
ADS
 
MATH
 
Google Scholar

Salomon-Ferrer, R., Gotz, A. W., Poole, D., Le Grand, S. & Walker, R. C. Routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit solvent particle mesh ewald. J. Chem. Theory Comput. 9, 3878–3888 (2013).

Article
 
CAS
 
PubMed
 
Google Scholar

Amber 16. University of California, San Francisco, CA, (2016).

Cornell, W. D. et al. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am. Chem. Soc. 117, 5179–5197 (1995).

Article
 
CAS
 
MATH
 
Google Scholar

Wang, J., Cieplak, P. & Kollman, P. A. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J. Comput. Chem. 21, 1049–1074 (2000).

Article
 
CAS
 
Google Scholar

Maier, J. A. et al. ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 11, 3696–3713 (2015).

Article
 
CAS
 
PubMed
 
PubMed Central
 
MATH
 
Google Scholar

Kirschner, K. N. et al. GLYCAM06: a generalizable biomolecular force field. Carbohydrates. J. Comput. Chem. 29, 622–655 (2008).

Article
 
CAS
 
PubMed
 
PubMed Central
 
MATH
 
Google Scholar

Lay, W. K., Miller, M. S. & Elcock, A. H. Optimizing solute-solute interactions in the GLYCAM06 and CHARMM36 carbohydrate force fields using osmotic pressure measurements. J. Chem. Theory Comput. 12, 1401–1407 (2016).

Article
 
CAS
 
PubMed
 
PubMed Central
 
MATH
 
Google Scholar

Lay, W. K., Miller, M. S. & Elcock, A. H. Reparameterization of solute-solute interactions for amino acid-sugar systems using isopiestic osmotic pressure molecular dynamics simulations. J. Chem. Theory Comput. 13, 1874–1882 (2017).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Miller, B. R. et al. (ed 3rd) MMpbsa.py: an efficient program for end-state free energy calculations. J. Chem. Theory Comput. 8 3314–3321 (2012).

Article
 
CAS
 
PubMed
 
MATH
 
Google Scholar

Roe, D. R. & Cheatham, T. E. III. PTRAJ and CPPTRAJ: Software for processing and analysis of molecular dynamics trajectory data. J. Chem. Theory Comput. 9, 3084–3095 (2013).

Article
 
CAS
 
PubMed
 
MATH
 
Google Scholar

Humphrey, W., Dalke, A. & Schulten, K. VMD: visual molecular dynamics. J. Mol. Graph 14, 33–38 (1996).

Article
 
CAS
 
PubMed
 
Google Scholar

The PyMOL molecular graphics system. ver. 2.5.0. https://pymol.org, New York, NY. (2021).

Download references

Acknowledgements

The quantum mechanics/molecular mechanics calculations were performed with funding from the Research Center for Computational Science, Okazaki, Japan (Projects 20-IMS-C027, 21-IMS-C024, and 22-IMS-C023), Promoted Projects by the Research Institute for Information Technology, Kyushu University, and the HPCI System Research Project (Projects hp220121, hp230171, and hp240062). We thank Dr Tim Cooper, from Edanz (https://jp.edanz.com/ac), for editing a draft of this manuscript.

Funding

This research was partly supported by JSPS KAKENHI (Grant No. JP21K05187 to T.Y. and Grant No. 23K18510 to T.U.) of the Japan Society for the Promotion of Science (JSPS).

Author information
Authors and Affiliations

Department of Applied Chemistry, Faculty of Engineering, University of Miyazaki, 1-1 Nishi, Gakuen-Kibanadai, Miyazaki, 889-2192, Japan

Toshifumi Yui & Takuya Uto

Contributions

T.Y. conceptualized the study, visualized the data, developed the methodology used, applied the software, performed the MD calculations, wrote the original draft, and reviewed and edited the manuscript. T.U. developed the force field parameters and developed the methodology used.

Corresponding author

Correspondence to Toshifumi Yui.

Ethics declarations
Competing interests

The authors declare no competing interests.

Additional information
Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2

Supplementary Material 3

Supplementary Material 4

Supplementary Material 5
Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article
Cite this article

Yui, T., Uto, T. Unbiased picture of the ligand docking process for the hevein protein–oligosaccharide complex. Sci Rep 15, 3335 (2025). https://doi.org/10.1038/s41598-025-87407-8

Download citation

Received
16 October 2024

Accepted
20 January 2025

Published
27 January 2025

Version of record
27 January 2025

DOI
https://doi.org/10.1038/s41598-025-87407-8

Share this article

Anyone you share the following link with will be able to read this content:

Get shareable link

Provided by the Springer Nature SharedIt content-sharing initiative

Keywords
Carbohydrate-binding protein
Chito-oligosaccharide
Unbiased ligand docking
Molecular dynamics simulation
Subjects
Computational biology and bioinformatics
Structural biology
Supplementary Material 2 Supplementary Material 3 Supplementary Material 4

Cancel

Preview · live