Dynamical Rearrangement of Human Epidermal Growth Factor Receptor 2 upon Antibody Binding: Effects on the Dimerization

Abstract

Human epidermal growth factor 2 (HER2) is a ligand-free tyrosine kinase receptor of the HER family that is overexpressed in some of the most aggressive tumours. Although it is known that HER2 dimerization involves a specific region of its extracellular domain, the so-called “dimerization arm”, the mechanism of dimerization inhibition remains uncertain. However, uncovering how antibody interactions lead to inhibition of HER2 dimerization is of key importance in understanding its role in tumour progression and therapy. Herein, we employed several computational modelling techniques for a molecular-level understanding of the interactions between HER and specific anti-HER2 antibodies, namely an antigen-binding (Fab) fragment (F0178) and a single-chain variable fragment from Trastuzumab (scFv). Specifically, we investigated the effects of antibody-HER2 interactions on the key residues of “dimerization arm” from molecular dynamics (MD) simulations of unbound HER (in a total of 1 µs), as well as ScFv:HER2 and F0178:HER2 complexes (for a total of 2.5 µs). A deep surface analysis of HER receptor revealed that the binding of specific anti-HER2 antibodies induced conformational changes both in the interfacial residues, which was expected, and in the ECDII (extracellular domain), in particular at the “dimerization arm”, which is critical in establishing protein–protein interface (PPI) interactions. Our results support and advance the knowledge on the already described trastuzumab effect on blocking HER2 dimerization through synergistic inhibition and/or steric hindrance. Furthermore, our approach offers a new strategy for fine-tuning target activity through allosteric ligands.

1. Introduction

The human epidermal growth factor receptor (EGFR) family comprises four members: human epidermal growth factor 1 (HER1), HER2, HER3, and HER4 [1]. They are overexpressed in a wide variety of malignant tumours, in particular in non-small cell lung cancer (NSCLC) (40%–80%), breast cancer (14%–91%), and head and neck cancer (80%–100%) [2]. In recent years, EGFR became an important biomarker and target of anti-tumour therapy [3]. Within the EGFR family, HER2 is one of the most important members and when overexpressed, represents one of the most aggressive phenotypes (approximately 20%–25% of invasive breast cancers) [4], making HER2 a well-established therapeutic target. Trastuzumab (Herceptin^®) was the first Food and Drug Administration (FDA)-approved antibody against HER2 and remains the gold standard for the treatment of HER2-overexpressing cancers [5]. Clinical studies showed that the combination of trastuzumab with adjuvant chemotherapy substantially improves disease-free survival (DFS) and overall survival [6]. From a structural point of view, HER2 is a 1255-amino-acid transmembrane glycoprotein with 185 kDa [7,8] known for its function as a ligand-less co-receptor. It is noteworthy that HER2 is the only member of the EGFR family that does not homodimerize under normal conditions, and its ectodomain does not directly bind ligands [9]. It is activated via heterodimerization with other ligand-bound receptors, with HER3 the preferred heterodimerization partner [10] or homodimerization when it is expressed at very high levels on the cell surface [11]. HER2 dimerization activates cytoplasmic kinase, which leads to autophosphorylation of its intracellular domains and thereby downstream signalling of cellular pathways. In normal cells, few HER2 molecules exist at the cell surface and, therefore, few heterodimers are formed, allowing growth signals to be relatively weak and controllable. However, if HER2 is overexpressed, the dimerization level increases, resulting in enhanced responsiveness to growth factors and malignant growth [8]. Structurally, HER2 comprises an extracellular domain (~600 amino acid residues, Figure 1A) with four extracellular domains (ECD) (ECDI–ECDIV), including two cysteine-rich domains (ECDII and ECDIV), a transmembrane domain, an intracellular tyrosine kinase domain, and a C-terminal tail with autophosphorylation sites [12].

The structures of the ligand-free HER2 ECDs (Figure 1A), ECDI and ECDIII, are very similar, with each containing a “β-helix” with sides, formed by three parallel β-sheets due to a high number of hydrophobic residues, mainly leucine. The spatial conformation of these two subdomains hinders the effective access in part explaining lack of ligand binding by HER2. ECDII and ECDIV also share a similar structure with small structural units held together by one or two disulphide bonds. ECDIV, a near transmembrane domain, is crucial to the stabilization of the protein–protein interactions between HER2 and its dimerization partner [15]. The structural conformation of ECDI and ECDIII brings them close together, allowing the ECDII projection to participate in receptor dimerization and initiation of signal transduction [16]. HER2 also includes a fixed, open conformation overhand outside of the ECDII to contact, in particular, with the so-called “dimerization arm” at Δ265–288 region of its dimerization partner. This region is a short hairpin loop [17,18] with conserved residues and plays an important role in the dimerization process [16,18,19]. Fu et al. have inferred the coupling mechanism of HER2 with the antigen-binding antibody fragment (Fab) F0178 based on the three-dimensional (3D) static crystal structure of F0178:HER2 [14]. The authors claimed that F0178 was able to bind to HER2 through a large and highly complementary interface including both ECDI and ECDIII expressing a high heterodimerization blocking activity [14]. However, the structural mechanism underlying HER2 dimerization is still unknown.

Here, we have employed a computational approach to deepen the current understanding of the binding mechanism of HER2 with specific anti-HER2 antibody fragments. For this purpose, we have selected two well-known anti-HER2-antibody fragments, namely F0178 and a single-chain variable fragment (scFv) from trastuzumab. Both antibodies inhibit HER2 with distinct mechanisms of action. The antibody fragment–HER2 complexes are depicted in Figure 1B. For the study of the ScFv:HER2 system we have considered the crystal structure of HER2 with trastuzumab Fab [13]. The crystal structure showed that the trastuzumab Fab binds to the C-terminal region of ECDIV of HER2, far away from the known dimerization region (ECDII). This interaction blocks proteolytic cleavage in HER2 and indirectly affects dimerization with other HERs [20]. It has been hypothesized that binding of trastuzumab could induce conformational changes providing a steric barrier avoiding kinase activation. This assembly could also promote steric blocking of HER2 by proteolytic cleavage of its extracellular region and thereby altering interaction between HER2 and other proteins [6].

In silico approaches are of crucial importance to clarify the structural and dynamical factors that influence protein–protein interfaces (PPIs). They allow the prediction of how these interactions occur, improving design and development of new therapeutic agents with better efficacies. Molecular dynamics (MD) approaches have been applied to the development and optimization of new probes with high affinity and specificity to HER2-positive tumours [21,22]. The calculation of binding free energies among mutants can be used to make a screening of the ligands library to be tested. Aiming at a detailed picture of the overall coupling mechanism between the different antibody fragments and the receptor as well as to understand the HER2 dimerization process, we employed computational modelling and molecular dynamics (MD) simulations of HER2 receptor in the apo form and upon complexation with F0178 and scFv from trastuzumab. We performed an undocumented systematic analysis of structural, energetic, and dynamic characteristics, ranging from pairwise interaction formation to covariance analyses for these complexes. Moreover, we identified the most important changes occurring at these PPIs and their influence on the conformational rearrangement of the specific domains involved in dimerization, which is key to the development of new HER2-specific therapies.

2. Materials and Methods

2.1. Model Construction

The 3D structure of HER2 extracellular domain was constructed by homology modelling using the MODELLER package [23], and the protein template NP_004439.2 (PDBID:3N85A) [24] with the sequence 301598662. The loop regions LEU99:SER133; THR245:CYS264; and ASP121:GLY136 were further refined with ModLoop [25,26]. The complexes of HER2 with antibody fragments were retrieved from the correspondent PDB files, specifically F0178:HER2 from PDBID:3WSQ [14] and ScFv:HER2 from trastuzumab with PDBID:1N8Z [13]. The protonation states of residues in the physiological pH range were determined with the PROPKA methodology [27,28,29], which computes pKa values of ionizable residues by accounting for the effect of the protein environment. PKA2PQR [27] was applied to convert PDB files to PQR molecular structure format files by AMBER ff99 force field [30]. All histidine residues were considered deprotonated at physiological pH conditions save for HIS 349, which was assigned as positively charged.

2.2. Molecular Dynamics Simulations

MD simulations of HER2, F0178:HER2, and ScFv:HER2 were performed using GROMACS 2018.3 [31,32,33] and the Amber ff99SB-ILDN force field [34]. Each complex was solvated by a TIP3P water box with 12 Å buffering distance to the boundary, which was verified throughout the simulations. The simulation systems were kept neutral by adding the necessary counterions. Simulations were performed in the NPT (isothermal–isobaric) ensemble. Temperature coupling was performed using the v-rescale thermostat [35] at 300 K, with a coupling constant of 0.1 ps, while an isotropic Parrinello–Rahman barostat [36,37] was used to keep the pressure constant at 1 bar using a coupling constant of 2.0 ps and a compressibility of 4.5 × 10⁻⁵ bar⁻¹. Long range electrostatic interactions were computed using the particle mesh Ewald (PME) method [38,39] with a Fourier grid spacing of 0.16 nm and a cutoff of 1.0 nm for direct contributions. Lennard-Jones interactions were computed using a nonbonded neighbour pair list with a cutoff of 1.0 nm, enabling the use of the Verlet scheme [40]. Solute bonds were constrained using the Parallel LINear Constraint Solver, P-LINCS [41], and solvent molecules were constrained using SETTLE algorithm [42]. The initial system energy was first minimized for 50,000 steps, followed by two equilibrations of 250 ps. The HER2 system was then initialized for 10 ps using position restraints on all protein heavy atoms with a force constant of 1000 kJ nm⁻² mol⁻¹, while the remaining systems were initialized for 25 ps without restraints. Initially, four production runs were performed for each system, each 200 ns long. In addition, for the F0178:HER2 and ScFv:HER2 complexes, a longer production run was also performed (~500 ns long). However, when analysing the centre-of-mass (COM) distance between the partners in the HER2:ScFv system (Figure S2), we decided to extend the length of all production runs to 500 ns, in order to obtain better sampling of this motion.

The equilibrated regions for each system (Table S1, in the Supplementary Materials) was replicate/system dependent, and the decision was based on the root-mean-square deviation (RMSD) and COM distance variations. The RMSD calculations were performed using the Cα atoms either in all extracellular domains, or in the ECDI–ECDIII region, as shown in Figure 2. The COM distance plots (Figure 3 and Figures S1 and S2) were obtained by calculating the centre-of-mass of the HER2 and either scFv or F0178, using all atoms in the receptor and each of its partners.

2.3. Solvent-Accessible Surface Area Analysis

We calculated the solvent-accessible surface area (SASA) of each individual residue in HER2, ScFv:HER2 and F0178:HER2 using the equilibrated regions of each replicate. In the unbound HER2, for instance, this resulted in four average values of SASA per residue, which were then used to calculate the mean and standard error of the mean (SEM) using a jack-knife technique. We then converted these absolute SASA values into percentages by calculating the maximum value a fully exposed residue of each type would have. With these percentages, we excluded all residues whose solvent access was below 20%, as they were considered mainly internalized.

2.3.1. SASA Analysis of Interfacial Residues

The identification of residues in the interface between HER2 and each of its partners was performed by calculating the SASA of each residue of the receptor in the complex, with and without its partner (e.g., the residue coloured in yellow depicted in Figure S6). When comparing the SASA values of each complex in the presence and absence of the antibody, we were able to exclude all residues whose SASA was not altered, i.e., residues which are not interacting directly with the partner. Finally, we only reported residues whose difference in SASA was above 5% and whose SEM was lower than this difference.

2.3.2. SASA Analysis of Allosterically Modulated Residues

To identify allosterically modulated residues, we compared the SASA of all residues in each complex with the same residue in unbound HER2. This SASA difference identified which residues became more exposed or occluded by an indirect effect of the antibody (e.g., the residue coloured in green depicted in Figure S6). As in the previous protocol, residues whose difference in SASA between systems was below 5% and whose SEM was higher than the difference itself were excluded.

2.4. Principal Component Analysis

Principal component analysis (PCA) [43] is a useful statistical technique able to reveal the more important patterns in a dataset. The most relevant information from a dataset is extracted as eigenvalues or principal components (PCs). Internal correlation motions between residues of interface complexes were investigated by performing a cross-correlations analysis on the MD trajectory given a covariance matrix in Cartesian coordinate space by Equation (1).

$$c_{ij}= \langle \left(x_i-\langle x_i\rangle \right)\left(x_j-\langle x_j\rangle \right)\rangle ,$$

(1)

where $x_i$ and $x_j$ are coordinates of the ith and jth atoms in the system and < … > denotes trajectory average [44]. The matrix $c_{ij}$ symmetric can be diagnosed by using orthogonal coordinate transformation matrix Ƭ, transforming the matrix $c_{ij}$ into a diagonal matrix Λ of eigenvalues λi by Equation (2):

$$\mathsf{\Lambda }={\mathrm{Ƭ}}^{\mathrm{Ƭ} }{c}_{ij}\mathrm{Ƭ} ,$$

(2)

where each eigenvector of the matrix gives the direction, along which a concerted motion arises, and the eigenvalue gives the magnitude of fluctuations along this direction. The PCA calculations were performed using bio3d [45].

2.5. Evolutionary Conservation of Interfacial Residues

Sequence conservation was retrieved for all interfacial residues in HER2 using ConSurf [46], which is based on the Rate4Site algorithm [47]. Multiple Alignment using Fast Fourier Transform (MAFFT) [48,49] was used for multiple sequence alignment (MSA), using BLAST [50] on the UNIREF90 database [51]. Each MSA comprised at most 150 sequences, with homology values ranging from 35% to 95%.

3. Results

3.1. Reaching Equilibration in the MD Simulations

To attain equilibrium in MD simulations of large protein complexes is not a trivial task. We need to evaluate the more conventional structural properties, such as the root-mean-square deviation (RMSD), but we will also need to monitor larger conformational transition of the quaternary subunits that may appear at the slower timescales. The RMSD values were calculated using the Cα atoms of the receptor, either using all ECD regions or excluding ECDIV (Figure 2). The inclusion of the ECDIV region in the calculations greatly increased the overall RMSD value due to its high conformational flexibility as can be observed in unbound HER2 and in the F0178:HER2 complex (Figure 2A,B). However, in ScFv:HER2, the contribution of ECDIV is reduced, which was assigned to the strong interaction of the scFv with this region (Figure 1B). Overall, the ECDI–III region equilibrates relatively fast, while the ECDIV RMSD values indicate larger movements occurring at larger timescales but of more difficult equilibration.

To investigate the slower conformational transitions between HER2 and its partners, we have also calculated the centre-of-mass (COM) distance between them. This distance is almost invariant in F0178:HER2 (Figure S1), since F0178 is interacting directly with the ECDI–ECDIII region of HER2 (Figure 1B), inducing no significant domain movement. In the ScFv:HER2 system (Figure S2), we observed that the scFv, upon binding, induced an opening/closing movement of the ECDIV towards the ECDI–III domains. As described in the methods section, this movement, which varies in the hundreds of nanoseconds timescale, was the reason why we have extended the MD simulations to improve its sampling.

Brought together, these results showed that scFv position varies between two extremes, depending on receptor proximity (Figure 3). The maximum COM distance leads to an opening of the cleft between scFv and the HER2 ECDI–III domains (Figure 3A), an intermediate position, which is similar to the initial conformation (Figure 3B); and a minimum distance, which leads to cleft closure and a direct interaction of scFv with ECDI–III (Figure 3C). Next, we selected the equilibrated regions of the MD simulations based on the RMSD and COM distances data. In the HER2 control and the F0178:HER2 complex simulations, these equilibrated rather fast (30–50 ns), whereas in the ScFv:HER2 system, due to the large domain movements, some replicates took up to 150 ns to reach equilibrium. In this system, all replicates seemed to converge to a different COM distance, indicating that the use of multiple replicates was crucial to improve the sampling of this movement.

3.2. HER2 Conformational Reorganization

After the equilibration evaluation of the MD simulations, we identified the equilibrium conformational space and performed a deeper analysis of HER2 conformation reorganization after antibody binding. Briefly, it has been previously described that HER2-ECDII is the main domain involved in the dimerization process [19,52]. Therefore, we have used an integrated approach of multiple methods to evaluate how HER2 coupling to its different partners could alter this functional mechanism. In addition, we have performed cross-correlation analysis for both complexes (Figure 4).

For the unbound HER2 (Figure 4A) we found a mixed correlation spectrum (both positive and negative) for the 23–186 region (ECDI) and negative correlation values (excluding ECDIV) for the rest of the protein. Fu et. al. identified HER2 sub-region ASN498–ASP502 (ECDIII) as particularly altered due to F0178 binding even though it is positioned slightly apart [14]. Cross-correlational analysis for F0178:HER2 shows strong and weak, respectively, negative correlation when comparing the 23–186 region (ECDI, especially sub-region 23–100) with ECDII (187–363) and ECDIII (364–530) regions. The negative correlation of ECDIII with ECDI (1–186 region) and ECDIV (531–626 region) underscores the importance of ECDIII on the F0178:HER2 interaction. Another very interesting aspect is the positive correlation between ECDIII and ECDII (187–363 region).

Concerning ScFv:HER2 (Figure 4C), a strong positive correlation between the ECDIV tail (above 800) and ECDI and ECDII revealed the importance of ECDIV on the ScFv:HER2 interaction. The strong negative correlation of ECDII and ECDIII, as opposed to HER2 movement when unbounded (Figure 4A) suggested a conformational change with scFv binding at ECDIV. This behaviour is in line with that shown in Figure 3 in which the ScFv:HER2 complex varies its position between two extremes.

To extract the relevant information from the conformational rearrangement observed during the MD simulations, we performed a PCA analysis for the F0178:HER2 and ScFv:HER2 systems. These calculations were performed using all replicas and considering only the equilibrated regions of the MD simulations. Conformational clustering was also used to identify and further illustrate the most populated conformers sampled in the MD simulations. For both systems, over 78%–95% of the movement was described by the first three principal components, which allowed us to positively identify the most relevant conformational transitions (Figure 5).

The PCA analysis for the F0178: HER2 MD simulations excluded the ECDIV tail, since its large movement (Figure 2B) would mask any other movement by the other domains. Assuming this approach allowed us to get insights about conformational details in regions otherwise inaccessible. Figure 5A shows slightly distinct conformational states mainly dependent on the ECDII at the 187–363 region. Remarkably, the ECDII (187–363 region) residues had a significant contribution to the first two principal components (Figure 5B). We also calculated the conservation scores using the ConSurf web server [53]. Figure S3 shows a high conservation score for GLU352, VAL353, and ARG354 at F0178:HER2 complexes, which is in accordance with the higher contribution of ECDII region for PC1.

PCA analyses suggested two main conformational states for the scFv:HER2 MD simulations (Figure 5C). In particular, the contribution of the ECDII residues to PC1 and PC2 (Figure 5D) suggested that the dimerization region was indeed affected by binding. Figure S3 also supports this behaviour with a high rate of the conservation pattern at ECDII, namely at dimerization arm (GLU265-ARG288). Paradoxically, the role of the ECDII-HER2 residues in the interface–interaction network of both complex systems were not considered in the literature [13,14]. To further aid interpretation, we produced a trajectory that interpolates between the most dissimilar structures in the distribution along PC1 for both systems (Figure S4). The movement of the “dimerization arm” at ECDII is rather evident, as well as at ECDI. The role of this region in the HER2 coupling mechanism has not been previously described; although the blocking ability of HER2 heterodimerization is well recognized [14]. We propose that the HER2 ECDII conformational change is coupled with this allosteric regulation of F0178:HER2.

Concerning the scFv:HER2 complex system, we observed an approximation between ECDII and scFv compared to the initial 3D crystallographic structure (PDBID:1NZ8 [13]) (Figure S4). This is an important biological feature since it is well known that the ‘heart-shaped’ domain II open conformation is essential for dimerization [17].

3.3. The Impact of ScFv and F0178 on the Conformational Space of HER2

Based on the analysis of the conformational changes that occurred during the MD simulations of the systems, scFv and F0178 antibody fragments seem to induce conformational changes in HER2, which might alter its dimerization propensity. The most impactful changes are those that alter the receptor surface, in particular, in the regions known to be involved in dimerization. These surface changes are inevitably coupled with variations in the solvent exposure of certain residues. Therefore, we calculated the solvent-accessible surface area (SASA) of each HER2 residue in all systems to calculate significant differences between complexes and the unbound receptor. The residues with significant differences in SASA (difference >5% and standard error is smaller than the difference value) due to the presence of a partner are associated with a direct interaction and should correlate well with the interfacial region (Figure 6). As expected, the residues with the largest shift in SASA in F0178:HER2 are located in ECDI, ECDII, and ECDIII (Figure 1B and Figure 6A), whereas in ScFv:HER2 these residues are mainly located in ECDII and ECDIV (Figure 1B and Figure 6B). In the case of ScFV:HER2, the direct interaction of the ligand with ECDII is due to the cleft closure (Figure 3C), observed in several simulations (Figure S2). The large error bars in F0178:HER2 were attributed to insufficient sampling and/or an increase variability in the interfacial region of the complex.

We also identified several residues whose SASA changed between unbound HER2 and its complexes through an indirect, or allosteric, effect of the antibody fragment. This was attained by calculating the SASA differences between F0178:HER2/ScFv:HER2 and HER2 trajectories, ignoring the physical presence of the ligands. The residues identified with this procedure (Figure 7 and Figure S5) were different in both complexes, but located predominantly in ECDII, including the so-called “dimerization arm” region (TYR274, GLU280, and MET282 in F0178:HER2; and CYS268 in ScFv:HER2).

Our results indicate that F0178 and scFv form fairly distinct dynamic complexes with HER2. The binding of F0178 seems to provoke a marked allosteric effect on the “dimerization arm” residues at ECDII, which is the dimerization region of HER2 [16]. A similar allosteric effect could be observed for scFv, which was complemented with a direct interaction with the ECDII region. This interaction map obtained from our simulations may help explain the strong effect of scFv on HER2 homodimerization [6], opening new therapeutic avenues for targeting HER2 inhibition.

4. Conclusions

HER2 plays a key role in the HER family due to its specificity to interact with other HER receptors through a complex signalling network to regulate cell growth, differentiation, and survival. HER2 overexpression occurs in approximately 15%–30% of breast cancers and 10%–30% of gastric/gastroesophageal cancers and serves as a prognostic and predictive biomarker. The easy accessibility to the four HER2 extracellular domains (ECD) turns this receptor into a highly relevant target for the selective delivery of anti-tumour drugs as well as imaging agents. ECDI and ECDIII can form a binding site for potential HER2 ligands, whereas the cysteine-rich ECDII and ECDIV are known to be involved in the homodimerization and heterodimerization processes. Furthermore, ECDII, in particular the “dimerization arm” region (Δ265–288), is believed to be the main contributor to dimerization. Although crystal structures are known, the dynamics of the binding of different antibodies to this receptor has remained unclear. In this work, we have applied an integrated framework of molecular modelling tools and techniques as well as MD simulation analyses to address important structural factors that provide atomistic information of the intermolecular coupling dynamics between HER2 and two different HER2 inhibitors, namely F0178 Fab and scFv from trastuzumab. Both are well-studied anti-HER2 antibody fragments that act at different ECD domains of HER2. F0178 interacts with both ECDI and ECDIII domains, while scFv binds at ECDIV domain. Concerning the ScFv:HER2 complex, the highest contribution to the PPI stability occurs in the terminal part of ECDIV.

In addition to the detailed residue and pairwise interaction analyses of HER2-antibody complexes, an important contribution of our work is on the conformational changes that occur upon binding. Our simulations suggested that F0178 and scFv form distinct dynamic complexes with HER2. The binding of F0178 seems to induce a marked allosteric effect on the “dimerization arm”, whereas scFv binding leads to major conformational changes in HER2 structure (ECDI approaches ECDII and ECDIII). We observed a drastic movement of ECDII towards scFv, establishing a new possible interface capable of meaningful and important interactions with the antibody. This newly discovered conformational organization of HER2 upon scFv coupling could indeed explain its role on dimerization as suggested by Hudis et al. [6], and herein we give a pioneer view of this crucial dynamical process. The integrated analysis framework applied here has unveiled the role of F0178 and scFv from trastuzumab in HER2 dimerization. Our results provide a new perspective for future approaches to the selection of effective potential drugs, since the conformational changes triggered by the ligands can have positive effects downstream of the tumour cell receptor pathways.