Introduction

The field of nanobodies has continuously developed over the past decades, and their characteristics and applications have gradually emerged. It took about 25 years from the serendipitous discovery of nanobodies to the FDA approval of the first nanobody drug, Caplacizumab1,2,3. Retaining full binding–antigen capacity, intensive studies on nanobodies have recently led to a wave of the “third generation” of antibodies, with growing interest in their treatment of diseases. In fact, the annual number of studies in this field increases exponentially4,5. In general, the spatial conformation of a nanobody contains two functional parts: the one includes four framework regions (denoted as FRs 1–4), and the other has three complementary determining regions (CDRs 1–3), as shown in Fig. 1. The CDR3 domain is highly flexible and primarily responsible for the binding ability of the coresponding nanobody6. Aside from their excellent property of small size (2.5–3.0–4.0 nm)7, nanobodies also exhibit high solubility, high stability8, and penetration9. In addition, nanobodies are characterized by either high affinities or high specificities in recognizing antigens and epitopes10,11,12. These remarkable features constitute the bases for their expanding application in various fields such as molecular imaging13,14,15, nanobody–drug conjugates16, drug delivery systems17, and especially in cancer diagnostics and therapeutics. Owing to their many promising advantages in both academic studies and commercial applications18, nanobodies continue to attract significant attention from the scientific community and have become a cornerstone in the field of immune-informatics19.

In the context of events occurring on the nanosecond time scale and at the atomistic level, recent computational advances provide us with powerful tools to enhance our understanding of protein–protein interactions. These advances give us a closer look into, among many others, the working mechanisms of macromolecular systems. Many algorithms based on rigid-body models such as the SDOCK by Zhang20and the MEGADOCK by Matsuzaki21, have been introduced to investigate protein–protein binding modes. Furthermore, the flexibility of both protein side-chains and main-chains has widely been incorporated into protein–protein interaction studies22,23. Moreover, the fast Fourier transform algorithm (FFT) used in ZDOCK24can exhaustively generate large numbers of decoys by rotating and translating the initial conformation of a given nanobody–protein complex. These developments have greatly inspired workers in the field to uncover nanobody–protein recognition25,26. Hacisuleyman27also claimed that the interaction energy is a distinguishing feature of protein–protein binding modes. The Dreiding interaction force field28(DI) energy method, when combined with ZDOCK, has proven quite useful. For example27, up to 35 out of 36 native forms of nanobodies have accurately been predicted using this theoretical approach. In 201729, the binding free energy of the EGFR inactive conformation with the nanobody 7D12 was estimated using the umbrella sampling method. Subsequently, Shao elucidated in 2019 the mechanical role of EGF in activating the transformation of EGFR in going from an inactive to an active state30. Computational approaches have proven to be effective not only in cancer research but also in addressing global issues such as the SARS-CoV-231,32. These advancements have greatly inspired us in approaching the structures and activities of nanobodies using theoretical methods.

The Epidermal Growth Factor Receptor, commonly denoted as EGFR or HER1, belongs to a family of receptors that also includes HER2, HER3 and HER433. EGFR is generally a part of all cells and is critical to the cell proliferation and embryogenesis34,35,36. EGFR is linked to cancer, as extensive evidence showed an overexpression of EGFR in cancer cells, namely more than 106EGFR molecules are found in cancer cells, whereas a normal cell contains approximately 30,000 receptors per cell37. Because this transmembrane protein is actually responsible for various types of tumors, it has become a primary target38in an anticancer therapy39. A typical wild-type EGFR contains three regions: an extracellular domain (containing 621 amino acids), a transmembrane domain (23 amino acids), and an intracellular region (542 amino acids). Two states of the EGFR extracellular domain have been identified, namely, the one has an active conformation and the other an inactive conformation. A cell signaling pathway begins with a transition of EGFR from its inactive to active state. Active EGFR monomers connect with each other to form homodimers (such as HER1 connects to HER1), or heterodimers (when HER1 connects to HER2, HER3, or HER4)40. To interrupt the downstream signaling pathway of EGFR, two basic strategies were identified: intracellular and extracellular inhibition. Indeed, the use of tyrosine kinase inhibitors (TKIs) for intracellular pathways and monoclonal antibodies for extracellular pathways has been proven to be effective approaches, already saving many human lives. Recently, due to the limitations of small-molecule TKIs41and antibody applications42, nanobodies have been introduced to enhance the treatment efficacy based on several advantages mentioned above including their small size, high penetration, specificity and selectivity43,44. Accordingly, the understanding on how a specific nanobody interacts with an epitope turns out to be one of the most important aspects of the current cancer research45.

To date, only a few structures of nanobodies46, including but not limited to the 7C1247,48,49, 7D1250,51, EgA152, EG253,54, EgC951, 9G855, 38G756, D10/E10/G1057and EgB458 have been reported in a bound state with the EGFR extracellular domain. Of the latter, only the EgB4 can target the active conformation of EGFR. The binding structure of the EgB4–EGFR dimer was isolated in 2022 by Zeronian et al.58 Although EgB4 has a shorter CDR3 region as compared to others (EgA1, 9G8, 7D12), its binding does neither compete with EGF nor induce any conformational change in the dimeric EGFR. The foreseeing application of EgB4 in development of drug conjugates or drug delivery systems necessitates an accurate evaluation of its binding affinity to the EGFR extracellular domain. In this context, we set out to explore the nature of interactions, and compute the binding affinity between the EgB4 antibody and the EGFR receptor that remain not well known yet.

The EGFR extracellular domain is a 170 kDa protein which can be split into four sub-domains including I (residues 1–164), II (residues 165–309), III (residues 310–479), and IV (residues 446–620)59. As for a convention, in what follows, residues of a nanobody are noted in bold lower case (such as asp98), whereas residues of EGFR protein are given in bold upper case (such as ARG141). Previous crystal structure analysis showed that the EgB4–EGFR complex is stabilized not only by hydrogen bonds between arg105 of the nanobody and LYS188, ILE189 and CYS191 residues at the top of EGFR but also by salt bridges between the nanobody asp98 and asp110 and the side chain ARG141of EGFR58. Moreover, amino acids such as trp53 and trp100 of EgB4, and THR140 and PHE156of EGFR collaboratively form a hydrophobic core. Intriguingly, the EGFR–EgB4 binding conformation is maintained, irrespective of the presence of EGF. The EGF ligand is found to play an important role in both homo- and hetero-dimerization of the EGFR extracellular domain60. While the binding site of EgB4 is located in domains I and II of EGFR, the binding region of the small EGF ligand is far from that, encompassing domains I and III. The EGF ligand interacts with the EGFR protein through both hydrophobic and electrostatic interactions via residues Y44, L47, and E5, D11, D17, R41, respectively30. Therefore, another effort in our present study is to detect the influence of the EGF molecule on the binding capacity of the EgB4 nanobody. This aspect has not well been understood yet.

In this study, even starting with an available crystal structure (7MO4), computational performance remains a quite challenging task. Performing massive MD simulations of a full-length EGFR dimer is rather difficult due to the large number of atoms in the system, in view of our actual limited computing resources. In a significant effort to approach the goals of this study, we prepare two solvated complexes containing a full-length monomeric EGFR extracellular domain and the EgB4 nanobody58; the one with the presence and the other without the presence of the small EGF ligand. Each simulated system thus contains up to a number close to 500,000 atoms.

First, we perform a set of classical molecular dynamics simulations to explore the EgB4–EGFR interaction. Several parallel trajectories are generated as initial conformations for each system to enhance conformation sampling. The most frequently appearing conformation, known as the representative structure, is extracted for further steps using steered molecular dynamics simulations and umbrella sampling method. Both the fast method (steered molecular dynamics) and the more accurate method (umbrella sampling) are applied to evaluate the binding affinity between EgB4 and EGFR molecules. Finally, our results predict that EgB4 binds to EGFR with more difficulty in the presence of EGF. To our knowledge, this is the first time the effect of EGF on the nanobody–EGFR binding is explored in detail.

Materials and methods

We use the available crystal structure 7OM4 which can be downloaded from the PDB Bank to prepare initial complexes between the EgB4 nanobody and EGFR protein, with one complex including the EGF ligand and the other without it58.

Molecular dynamics simulation methods (MD)

All molecular atoms are generated using a Gromacs version 2020’s gmx pdb2gmx module61. The water model TIP3P and the force field amber99sb-ildn are selected62,63. Before adding water molecules and charge-neutralizing ions, the entire complex structure is arranged in a cubic box in such a way that the distances between the complex and the box borders are greater than 12 angstroms. This results in each system containing a number close to 500,000 atoms. After that, the system reached an equilibrium in three stages, an energy minimization using the steepest descent algorithm64, a 10 ns constant volume (NVT), and 10 ns constant pressure-temperature (NPT) equilibrium. The P-LINCS algorithm is used to handle all hydrogen bond lengths65. For both the Van der Waals and electrostatic interactions, a cutoff radius of 10 Å is used. Electrostatic interactions are calculated using the particle mesh Ewald method66. Pressure is maintained at 1 bar using Parrinello-Rahman67 pressure coupling, and temperature is kept at 298 K using velocity rescaling with a stochastic term. Either the EGF ligand is present or absent, five independent trajectories and 500 ns of simulation time per trajectory are carried out to enhance conformation sampling. Representative structures extracted from MD steps are used as inputs for both SMD and PMF calculations.

Steered molecular dynamics simulation

In a typical steered molecular dynamics (SMD) simulation, an external force is applied to the center of mass (COM) of EgB4, driving the nanobody away from its initial position68,69,70. The pulling velocity and spring constant are set to 5 m/s and 600 kJ/mol/nm², respectively. The pulling direction is shown in Fig. 1. The simulation time is set to 1 ns for all trajectories. In each pulling trajectory, snapshots are recorded every 2,000 fs, and the force–time and displacement–time profiles are collected every 10 fs. Notably, to ensure the external force breaks the EgB4–EGFR association, we restrain the EGFR backbone with a weak harmonic potential of 1,000 kJ/mol/nm². In this study, all Ca atoms of domains III and IV are fixed to keep the EGFR protein more flexible. A total of 200 pulling trajectories are performed for each system, the one with the presence of EGF and the other without it.

Umbrella sampling method

Before running the umbrella sampling method, steered molecular dynamics (SMD) simulations are employed to collect a group of EgB4–EGFR configurations because the US method, by its characteristics, cannot generate a dissociation process. In this pulling simulation, all EGFR Ca atoms are restrained, the pulling velocity is set to 1 m/s to slow down the escape of the nanobody, and the simulation time is increased to 4 ns. During the unbinding process, more than 70 windows are collected on the basis of the movement of the EgB4 center of mass by 0.05 nm from adjacent windows. Each of these 70 windows is then subjected to a 10 ns MD simulation. Due to the large size of the complex, containing nearly 500,000 atoms, we could only replicate two independent works for each condition, with both the presence and absence of the EGF ligand. The outcome of the umbrella sampling is finalized using the WHAM package of Gromacs 2020 to extract the PMF curves.

Hydrogen bonds and contacts

Formation of hydrogen bonds between any potential donor and acceptor is examined using the gmx hbond module in the GROMACS package. A hydrogen bond is considered to occur when the following geometric criteria are satisfied: (1) the distance between the donor (D) and acceptor (A) is less than or equal to 3.5 Å, and (2) the donor–acceptor angle is less than or equal to 30°. Note that the distance of 3.5 Å corresponds to the minimum of the radial distribution function (RDF) of water molecules in the SPC water model. On the other hand, a contact is formed in this study when an atom of the nanobody and an atom of EGFR are within 0.6 nm of each other. The number of contacts is counted using the gmx mindist module.

Interaction energy

The non-bonded interaction energy between two objects is often the sum of two terms including (1) the van der Waals (vdW) interaction energy provided by the Lennard-Jones (LJ) potential and (2) the electrostatic interaction energy provided by the Coulombic potential. To calculate this quantity, we need to split the full set of MD trajectories into several groups of snapshots, and then run a one-step molecular dynamics simulation on each snapshot. In the present study, the energy groups between the nanobody and EGFR are then defined.

Fig. 1
figure 1

The overall structure of EGFR (in green – cygan – wheat and orange), EGF ligand (in red), and EgB4 nanobody (in magenta). Zoomed-in view of the interface between domain I and II (from PRO130 to CYS208) of the EGFR protein (in green – cygan) and the CDR1 (in yellow), CDR2 (lightorange), and CDR3 (in blue) regions of the EgB4 nanobody. The pink cartoon represents the frameworks of the nanobody.

Full size image

Results and discussion

Interaction dynamics between EgB4 and EGFR over time

Figure 2 shows the root mean square deviations (RMSDs) obtained from our 10 µs molecular dynamics simulation over time. The RMSD values for all atoms of EgB4 averaged between 0.04 nm and 0.2 nm are listed in Table 1. The RMSDs of EgB4 remain consistent across all five trajectories, consistent with the high stability characteristic widely known for nanobodies.

Fig. 2
figure 2

Root mean square deviations (RMSDs) of all Ca atoms of EGFR and EgB4 nanobody in the presence and absence of the EGF ligand.

Full size image

When we examine the RMSD of EGFR for all atoms, in the presence of the EGF ligand, the mean values averaged from the first 300 ns and the last 200 ns of five trajectories amount to 0.11 nm and 0.12 nm, respectively. Conversely, in the absence of the EGF ligand, the RMSD values for the EGFR protein become nearly three times as large, measuring 0.59 nm and 0.65 nm for the first 300 ns and last 200 ns of the five MD trajectories, respectively.

During each trajectory having a 500 ns duration, the time-dependent number of hydrogen bonds and the interaction energy between the EGFR protein and the EgB4 nanobody are calculated. The EGFR-EgB4 complex formed an average of 6.32 hydrogen bonds per frame during the first 300 ns and 6.25 hydrogen bonds during the last 200 ns in the presence of EGF, as listed in Table 1. In contrast, when EGF is absent in the interacting system, these values are getting notably larger, averaging as 7.23 and 7.73 hydrogen bonds per frame during the first 300 ns and last 200 ns, respectively. This illustrates that most of the time, in the absence of EGF, formation of the EGFR-EgB4 complex induces more hydrogen bonds than when EGF is present.

In addition to the hydrogen bond profiles, the interaction energies between the EGFR protein and the EgB4 nanobody plotted over time are shown in Fig. 3.

Table 1 Numerical data including RMSD, hydrogen bonds, electrostatic and Van Der Waals potentials and interaction energy, averaged from five independent trajectories of EGFR–Egb4 in the presence and absence of the EGF ligand. Data are separated into two time periods: the first 300 ns and the last 200 ns of molecular dynamics simulations.
Full size table

No significant change is detected between the averaged values of interaction energy from the first 300 ns and the last 200 ns of each 500 ns EGFR–EgB4 simulation. In the first 300 ns, when EGF binds to EGFR, the interaction energy averages to −139.7 kcal/mol, while in the last 200 ns, it converges to −140.0 kcal/mol. Conversely, in the absence of EGF, these values become substantially larger, averaging at −151.6 kcal/mol and − 157.1 kcal/mol.

This numerical MD result highlights the role of the EGF ligand. The presence of EGF thus tends to enhance the flexibility of the EGFR protein. Moreover, the interaction energy between EGFR and the EgB4 nanobody is primarily caused by electrostatic interactions, increased from − 80.8 kcal/mol to −92.1 kcal/mol with the presence of EGF, and from − 80.6 kcal/mol to −95.4 kcal/mol without EGF. Further details regarding the contribution of each residue are described in the following paragraph.

Fig. 3
figure 3

The averaged time-dependent number of hydrogen bonds per frame, interaction energy, Coulomb potential, and van der Waals potential between the EGFR protein and the nanobody EgB4 were obtained from 5 independent trajectories of 500 ns MD simulations in the presence (green) and absence (red) of EGF molecules.

Full size image

Diversity of hydrogen bonds created by key residues of EGFR and EgB4

Twenty-two important residues that constitute the CDR2 and CDR3 regions of the EgB4 nanobody, are examined in this section. Particularly, during the last 200 ns of all simulations, the EGFR – asp98 and EGFR – asp110 salt bridges form an average of 1.2 and 0.9 hydrogen bonds per frame in the presence of the EGF ligand, and 1.0 and 1.5 hydrogen bonds per frame in its absence. An identical set of values is obtained when we count the number of hydrogen bonds with the chosen options, only ARG141 and asp98 / asp110 (data not shown). Hence, it is recognized that both asp98 and asp110 of EgB4 do not engage in a donor-acceptor relationship with any residues of EGFR, except for ARG141. In contrast to the reduction of hydrogen bonds involving asp98, the residue asp110 of EgB4 increases its hydrogen bond formation with ARG141 when EGF is absent.

A similar observation is made for hydrogen bonding between asp108 and ARG141. The involvement of asp108 in the association of EgB4 and EGFR has not been reported in previous experimental studies. Actually, asp108 of EgB4 is a highly active residue as it increases the average number of hydrogen bonds per frame with ARG141 from 0.14 to 0.64 in the presence and absence of EGF, respectively. The activity of asp108 is difficult to be detected without computational support.

Table 2 shows the average values obtained from five independent trajectories. No hydrogen bonds are found between seven residues including thr50, ile51, ala97, ala101, ser103, tyr109 and tyr111 with the EGFR receptor. Nine residues of EgB4 formed either a small number of hydrogen bonds (ser52, thr54, arg99, arg104) or an extremely small number of hydrogen bonds (thr53, asp55, ser56, asn106, val107) irrespective of the presence or absence of the EGF ligand. According to the results plotted in Fig. 4, when the system is changed from the presence to the absence of EGF, an increase is observed in the number of hydrogen bonds for asp108 and asp110, and a decrease for asp98, trp100, ser102, and arg105. During the last 200 ns of the five trajectories, the average number of hydrogen bonds formed between EGFR and trp100, ser102, and arg105 are equal to 0.65, 0.73 and 0.69, respectively.

Fig. 4
figure 4

Average number of hydrogen bonds per frame formed between 6 key residues of EgB4 and the EGFR protein. Data were obtained from the last 200 ns of five independent trajectories in the presence (in green) and absence (in red) of the EGF molecule.

Full size image

Information about how these amino acids interact with EGFR is crucial for rationalizing the stabilization of the EgB4-EGFR association. Besides ARG141 which is highly interactive, SER137, GLN139, ASN172 and GLN193 are found among the most active residues of the EGFR protein, frequently leading to formation of various bonds with either CDR2 or CDR3 of EgB4. Table S1 (SI file) points out the fact that the absence of EGF allows more residues of the EGFR receptor to contribute to the diversity of hydrogen bonding.

Table 2 Average number of hydrogen bonds, electrostatic potential energy and VdW potential energy per frame formed by key EgB4 residues and the entire EGFR protein. Results are estimated from the last 200 ns of five independent trajectories in both the presence and absence of the EGF ligand.
Full size table

Interaction energy

To better understand the role of each residue in the interaction energy between the nanobody and EGFR, the Coulombic and Van der Waals potentials between key residues of EgB4 and EGFR are evaluated. As detailed in Table 2, compared to others, five important residues including asp98, trp100, arg105, asp108 and asp110 are characterized by the greatest interaction energies of −24.0, −19.2, −19.8, −10.8 and − 21.9 kcal/mol, respectively, when EGF is present, and − 19.3, −17.9, −15.3, −15.5 and − 28.3 kcal/mol, respectively, when EGF is absent. These values explain for why these residues generate larger numbers of hydrogen bonds, as mentioned in the previous paragraph. Particularly, the electrostatic potentials of asp98 and asp110 contribute dominantly to the interaction energy, with values of −22.2 and − 21.1 kcal/mol in the presence of EGF and − 18.7 and − 22.6 kcal/mol in its absence. These values are much larger than the corresponding Van der Waals potentials of −1.4, −0.8, −0.6, and − 1.7 kcal/mol, respectively.

On the other hand, trp100 strongly interacts with EGFR via the Van der Waals force, with averaged values of −14.6 and − 14.8 kcal/mol collected from the last 200 ns of the two systems with and without EGF. For trp100, electrostatic interaction is relatively small, at −4.6 and − 3.1 kcal/mol, respectively. In the case of asp108, the electrostatic energy increases from − 6.2 to −12.2 kcal/mol, causing a difference of −4.7 kcal/mol in interaction energy, despite its Van der Waals energy decreases from − 4.4 to −2.9 kcal/mol. The largest interaction energy difference of −6.4 kcal/mol belongs to asp110.

In summary, while the absence of EGF enhances the interaction energy between asp108, asp110 and EGFR, it restrains the interaction of asp98, trp100, and arg105, resulting in positive potential differences of 4.7, 1.3, and 4.5 kcal/mol, respectively (cf. Table 5). Interaction energies from other residues converge to smaller average values for ser52, thr53, thr54, ser56, arg99, ala101, ser102 and arg104 or even smaller than − 3.0 kcal/mol such as the case of thr50, ile51, asp55, ala97, ser103, asn106, val107, tyr109, and tyr111.

From a side-view of EGFR, we also compute the polar/non-polar interaction energies between 97 EGFR residues located in the interface and EgB4 nanobody Fig. 5. Full data are listed in Table S3 (SI file). Besides the contributions from GLN139, ASN172, ILE189, our data analysis proves the key role of ARG141 when significant values of total interaction energies either in the presence (−57.5 kcal/mol) or in the absence (−57.2 kcal/mol) of EGF are computed.

Fig. 5
figure 5

Electrostatic energy, Van der Waals energy and interaction energy between EgB4 and EGFR in both the presence and absence of the EGF ligand.

Full size image

Representative conformations of EGFR-EgB4 complex in the presence and absence of the EGF ligand

In this section, the total number of hydrogen bonds and total interaction energy following interaction between both EgB4 and EGFR entities are used as variables to extract their representative structures. Figure 6 displays the snapshot corresponding to the location with the highest distribution (denoted as the global minimum) which is chosen as the representative conformation. Given that a snapshot is saved every 5,000 steps of the MD simulation, data from the last 200 ns of each trajectory generated a total of 10,000 snapshots for each case of EGF presence and absence. As a result, Fig. 6 shows the free energy landscapes and interaction pictures in both situations. The most important interactions are indicated. The structural comparison between representative conformations and initial conformation is also showed in the Figure S15 (SI file). The representative structures obtained in this part are subsequently used as initial structures for pulling simulations and umbrella sampling methods.

Fig. 6
figure 6figure 6

Free energy landscape of EgB4-EGFR systems in the presence of EGF (upper panel) and absence of EGF (lower panel). Representative structure is also shown in lower panel. Hydrogen bonds are indicated by green dotted lines and salt bridges by red dotted lines.

Full size image

Absence of the EGF ligand expands the interaction space between EgB4 and EGFR

In SMD simulations, to ensure that the external force breaks the EgB4-EGFR association and to elucidate the role of the EGF small ligand, we maintain the flexibility of the EGFR backbone. Therefore, only a small group of Ca atoms is restrained as described above in the Methods section. This approach, previously employed in the context of proteins, treats them as dynamic breathing objects and ensures that the applied force not only focuses on breaking the EGFR-EgB4 association but also induces additional work, potentially stretching the EGFR backbone whether the EGF is present or not. The representative root mean square deviations of the EGFR backbone are illustrated in Figure S6 (SI file).

Fig. 7
figure 7

Number of contacts between asp98 and ARG141 formed as a function of displacement of EgB4’s center of mass. Similar images can be viewed in Figure S5 (SI) file.

Full size image

Let us note that a comparison between the rupture force and external work is not reasonable, because the results obtained are not consistent enough to draw any conclusive consensus. Moreover, an accurate estimation of this part of the pulling work is not fully resolved yet. Therefore, a certain relationship between these typical non-equilibrium quantities cannot further be discussed. The time-dependent forces and works are presented in Figures S1-S6 (SI file).

Table 3 Average distance the EgB4’s center of mass moves to break contact with specific residues in the presence and absence of the EGF ligand. Data are obtained from 200 independent trajectories in each case.
Full size table

Another aspect of interest focuses on the calculation of the average number of contacts per frame, based on the distance that the EgB4 nanobody is moved. Numerical results obtained reflect the breadth of the EgB4-EGFR potential energy landscape. From previous sections, it is evident that seven key-pairwise interactions predominantly contribute to the EgB4-EGFR association, including the asp98-ARG141, trp100-ARG141, trp100-TRP140, ser102-GLN139, arg105-CYS191, asp108-ARG141 and asp110-ARG141. In referring to Fig. 7 and S5 (SI file), it is apparent that in the absence of EGF, all seven pairs maintain contact with each other over greater distances as compared to the situation when EGF is present. Averaged results taken from 200 independent trajectories are detailed in Table 3.

As the nanobody moves, with an EGF bound, under an external pulling force away from EGFR interaction, each pair of contacts disappears after distances ranging from at least 0.55 nm (ser102 – GLN139) to up to 1.12 nm (trp100 – ARG141). Conversely, in the absence of EGF, the smallest and largest distances observed before disruption of contacts are equal to around 0.77 nm (arg105 – CYS191) and 1.48 nm (trp100 – ARG141), respectively. This intriguing observation suggests that the absence of EGF tends to increase the volume of host space within which both EgB4 and EGFR entities can interact with each other.

Binding ability of EgB4 to EGFR

Free energy variation (ΔG) is a critical parameter in the study of biological systems. We now construct the free energy profiles, as shown in Fig. 8, through extensive efforts. In total, over 280 windows are used in this approach with two diagrams shown in Figure S8 (SI file). The convergence of umbrella sampling method is evaluated and presented in the Figures S9, S10 and S11 (SI file). Along the reaction coordinate, as the EGFR-EgB4 complex system undergoes a transition from a bound state to an unbound state, free energy values of −13.3 and − 17.1 kcal/mol are obtained with and without EGF binding, respectively. In reality, an enhancement in binding free energy occurs when the EGF is absent.

Fig. 8
figure 8

Equilibrium free energy profile of the EgB4-EGFR system in the presence (green) and absence (red) of the EGF ligand. Dashed arrows indicate the free energy difference between the bound state and the unbound state.

Full size image

Concluding remarks

In Figure S12 (SI file), the absence of EGF induces an initial space between both domains I and III. Logically, this tends to drive the system into two different situations, namely, either decrease the space or increase the space between 2 domains. Figure S13 (SI file) presents the number of contacts between two domains with respect to the evolution of time. It shows that:

  1. i)

    two domains I and III tend to attract to each other in the case of EGF absence and.

  2. ii)

    none of EGFR conformation is changed in the case of EGF presence.

In reality, the “active” EGFR protein without an EGF ligand has been known that it is easily transformed into an “inactive” state (showed in Figure S14, SI file). However with 500ns MD simulations of each trajectory, it is apparent that our simulations are not long enough to achieve the final motions of 4 domains. In addition, the influence of the nanobody EgB4 into the EGFR structural transformation remains not clear.

All these questions need to be resolved not only by extending the simulation work into micro-second time frame (even longer), but also by simulating more systems such as the EGFR protein without EgB4 nanobody. It will in fact consume much bigger computing resources that go beyond our current capacities. Therefore, in the scope of the present study we cannot take any definitive conclusion about the EGFR structural change, irrespective of the presence or absence of the EGF ligand. We would suggest considering this aspect in following studies.

Inspiringly, although the results related to the motion of domains are not fully adequate, we find that the interface between nanobody and EGFR is stable after 5 trajectories of 500ns MD simulations when we construct the free energy surface in Fig. 6 from the number of contact (PC1) and the interaction energy (PC2) between EgB4 and EGFR. Representative conformations of the EGFR-EgB4 interface are displayed in Figure S15 (SI file).

Across a total of five microseconds of simulations involving interacting systems containing up to 500,000 atoms, we extensively study the interactive dynamics between the EGFR protein and the EgB4 nanobody at an atomistic level. Significantly, numerous hydrogen bonds and interactions are identified both in the presence and absence of an EGF ligand. The contribution of asp108 to the association between EgB4 and EGFR protein is notable, alongside with the asp98, trp100, ser102, arg105 and asp110. In the absence of the EGF ligand, the role of asp108 is markedly increased. EGFR residues including GLN139, ASN172, ILE189 and especially ARG141 also hold the main role in the interaction dynamic with EgB4.

Our study also aims to investigate on how the EGF ligand affects the binding ability of EgB4 to the EGFR receptor. Calculated RMSD profiles from both MD and SMD trajectory runs indicate that the absence of an EGF ligand induces greater fluctuations in the EGFR structure. Our findings also suggest that the absence of EGF not only allows the EGFR to fluctuate more freely but also enhances its binding affinity to the EgB4 nanobody. Unfortunately, according to previous reports, approximately 95% of EGFR dimeric conformations are found in a binding mode with the EGF ligand71,72. In other words, when EGF binds to and stabilizes the dimeric structure of EGFR, lower opportunity is present for the nanobody EgB4 to bind to EGFR. This poses a significant challenge if the EgB4 nanobody is intended for use as a binder of EGFR in practical applications such as in biosensor development.

However, the present study also suggests another strategy for further exploration, that is, mutating the EgB4 structure based on the contact map plotted in Figure S7 (SI file). Specifically, we can consider a replacement of residues in the CDR3 region that either do not interact with EGFR or have very weak interactions, such as the arg99, asn106, val107 and trp109. By pursuing this strategy, a stronger anti-EGFR nanobody can be further developed.