Introduction

The global impact of the COVID-19 pandemic has led to extensive research on prevention of infection through development of effective vaccines and inhibitors. While some groups have taken an experimental approach to search for novel inhibitory molecules, many others have used computational methods due to their time and cost-effectiveness (De Vivo et al. 2016). The research on inhibitory molecules for SARS-CoV-2 has been directed to two main targets, the S1 domain of spike (S) protein on the virus and the angiotensin-converting enzyme 2 (ACE2) protein on the host cells, which serves as a binding site for S protein (Sawant et al. 2021; Srivastava et al. 2021; Ridgway et al. 2022). In computational studies, a wide array of potential inhibitors are docked to either S1 or ACE2 and are then evaluated based on the thermodynamics and/or the kinetics of their interaction. While computational methods can allow researchers to evaluate many potential inhibitors in a relatively short period of time, selection of appropriate computational models are critical to generate results that closely reflect the experimental data.

Many receptors and enzymes contain metals ions such as zinc and copper, serving an important role in the structure and/or function of these protein molecules (Banerjee et al. 2016; Putterill et al. 2022). It is estimated that nearly half of all human proteins contain at least one metal ion (Sakharov and Lim 2005). These metal ions are often essential for catalytic or co-catalytic activity, and can have large implications on the overall structural stability of the protein (Auld 2001; Zughaier and Cornelis 2018; Hübner and Haase 2021). Various functional groups within the protein can coordinate to metal ions, which form strong bonds that increase the structural stability within the region. Additionally, metal ions can have far-reaching effects on the movement of distant atoms within the protein due to interactions through hydrogen bonds (Christianson 1991). Despite their abundance and importance in protein structure, metal ions can be difficult to accurately parameterize for MD simulations. Metal-bound proteins can contain complex coordination geometries with strong protein-metal interactions, charge transfer, polarization, and even varying coordination numbers (Sakharov and Lim 2005). Although the use of bonded or non-bonded models to estimate the coordination geometry are quite uncommon, it may simplify the computational modeling of the given metalloprotein. Additionally, excluding the metal ions may also affect the overall protein structure which could cause computational results to deviate from the experimental observations.

The structure of ACE2 has some complexities that can make it difficult to accurately model. In particular, ACE2 contains a Zn ion at the angiotensin catalytic site, which is approximately 22 Å away from the S1 binding site. Zn is well-known for being difficult to model accurately due to its unusually strong and complex interactions with select amino acids (Laitaoja et al. 2013). While few previous studies have explored parameterization and computation of the Zn ion, due to both difficulties in modeling and its distance from the S1 binding site, many computational studies have chosen to exclude Zn from their models of ACE2 (Basu et al. 2020). However, studies have shown that Zn can have large effects as far away as 34 Å, so its exclusion may negatively impact the accuracy of these models (Sinha et al. 2016).

The computational structural immunology is an emerging field useful to predict effective viral inhibitors and drug targets. The conformational study of Zn-metalloprotein helps to better understand the structural changes due to Zn binding and to assess how they may impact proteins stability and function (Chaturvedi and Shrivastava 2005; Krężel and Maret 2016; Greener et al. 2018). In our previous work, we have explored the implications of Zn-bound angiotensin peptides on ACE2-S1 binding (Fatouros et al. 2022). Additionally, the time-based dynamics of immunologically relevant protein systems were discussed (Roy 2016, 2019). In this study, the authors evaluate the importance of Zn within model structures of ACE2-S1 and a monoclonal antibody (mAb) inhibitor-based ACE2-inhibitor complex by examining its impact on structural stability and binding free energies. As metalloproteins are important for ligand binding and various bio-catalytic processes, it is critical to identify essential non-protein features to ensure that computational models remain efficient yet accurate in drug design.

Materials and methods

Overview

To examine the differences attributed to the presence of Zn, crystal structures for ACE2-S1 and ACE2-inhibitor complexes were obtained and Zn ions were parameterized/excluded, followed by MD simulations for both Zn-bound and Zn-free versions of the complex. Structural stability and binding free energy were then computed and analyzed using the results of the MD simulation.

Modeling ACE2-S1 complexes

The experimentally derived structure for the ACE2-S1 complex, chains A and E from 6M0J (Lan et al. 2020), was obtained from RCSB PDB. Chain E is denoted in this report as S1. The Zn-free model structure was created by deleting Zn from one copy of the structure. Amber parameters, neutralizing ions (Na+ and Cl−), and a water box were generated for the Zn-free complex using the LEaP program (Case et al. 2020).

The structure and docking configuration of an ACE2 inhibitor, monoclonal antibody (mAb) chains H and L from 7E7E (Du et al. 2021), was also obtained from RCSB PDB. Although ACE2 is included in this experimentally derived structure, its structure is broken by a section of missing residues. Therefore, to replace this chain with an unbroken ACE2 structure, chains H and L of 7E7E were docked to chain A of 6M0J on the LZerD Protein Docking Webserver (Christoffer et al. 2021). To increase the chances of obtaining structures similar to the ACE2-inhibitor complex (henceforth designated as ACE2-mAb complex since the inhibitor is a mAb in nature) in 7E7E, all residues within 7 Å of the binding interface were specified in the docking program, although the allowed distance was increased to 10 Å to account for slight positional differences in ACE2. The top 500 docked structures were compared to the experimentally derived ACE2-inhibitor complex, 7E7E. The docked complex with the lowest root mean square deviation (RMSD) from the position of the experimentally derived structure was selected for analysis. Similar to the ACE2-S1 Zn-free complex, the LEaP program was used to add neutralizing ions, a water box, and parameterize the system.

In both systems, the Zn-bound protein complex was created by parameterizing Zn. To parameterize Zn with the appropriate coordination geometry, Amber force field parameters were generated for Zn and each coordinated residue using the bonded model in MCPB.py (Li and Merz 2016). In this method, force field parameters were generated with the Empirical method. During restrained electrostatic potential (RESP) charge fitting, the charges on the heavy backbone atoms of coordinated residues were constrained. Additionally, to convert from CHARMM to Amber notation and estimate residues’ protonation states at a pH of 7.0, Bio3D (Grant et al. 2006) and H++ (Gordon et al. 2005) were used sequentially. Density-functional theory was also employed, via GAMESS-US (Barca et al. 2020), to calculate the charges of atoms within the coordinating residues. Similar to the Zn-free complex, the LEaP program was used to add neutralizing salt and a water box with 10 Å before molecular dynamics simulations (Case et al. 2020).

Molecular dynamics simulations

After parameterizing both the Zn-bound and Zn-free complexes, MD simulations were performed in NAMD (Phillips et al. 2005) to minimize energy and then to explore structural differences. Minimization was performed using a step size of 1 fs at 310 K, controlled through Langevin dynamics with damping coefficient of 5 ps−1, and 1 atm, controlled with a Langevin piston with an oscillation period of 100 fs and damping timescale of 50 fs. To ensure that both total protein RMSD and total energy were achieved, VMD plugins RMSD and NAMD Energy were used (Humphrey et al. 1996). While the ACE2-S1 complex required 30,000 iterations for minimization, the ACE2-inhibitor complex required 60,000 iterations due to its greater size and instability from the added docking step. A constant volume MD simulation of each complex was performed with a step size of 1 fs at 310 K, controlled through Langevin dynamics. A constant volume equilibration MD simulation was followed by a 50 ns MD production run. For the smaller ACE2-S1 complexes, a 2 ns equilibration was required to stabilize RMSD, while a 10 ns equilibration was required for the larger ACE2-inhibitor complexes. Adequate equilibration was defined by stabilization in both protein RMSD and the system’s total energy. For all MD simulations Particle mesh Ewald (PME) was used to account for electrostatic interactions.

Analysis of structural stability and binding free energy

Once the 50 ns production runs were complete, they could be used to analyze the structural stability and binding free energy of both the Zn-bound and Zn-free complexes. To evaluate structural stability, RMSD, root mean square fluctuations (RMSF), and secondary structure of each complex during the 50 ns MD simulation were visualized with VMD plugins. The binding free energy of each protein complex was calculated employing the MM/PBSA method, using the CaFE plugin for VMD (Liu and Hou 2016). To calculate the solvent accessible surface area (SASA), APBS was used (Baker et al. 2001). Due to the system size and length of the MD simulation, binding free energy was computed with a step size of 0.45 ns. The solvent dielectric constant was kept at 80 while internal dielectric constants of 1, 2, 4, and 6 were used. As previously described by Fatouros et al., larger internal dielectric constants are required for polar or charged binding sites (Fatouros et al. 2022).

Results

The potential changes in the structural stability of the ACE2-S1 and ACE2-mAb complexes were explored by computation of RMSD and root mean structure fluctuation (RMSF). In addition to the overall complex, measurements were made for the interacting residues (IR), defined as the residues within 3.5 Å of the opposite protein within the unmodified 6M0J and 7E7E crystal structures. From Fig. 1A, it appears that the Zn-bound ACE2-S1 complex has greater average RMSD and RMSF than its Zn-free counterpart, while the IR exhibit lower RMSD in the Zn-bound complex. This suggests that the ACE2-S1 Zn-bound complex is less stable overall, but more stable at the IR. In contrast, within the ACE2-mAb system the Zn-bound complex has lower average RMSF and RMSD, both in the overall protein and at the IR (Fig. 1B). As anticipated, ACE2, which contains the zinc ion, shows a greater difference in average RMSF in both the ACE2-S1 and ACE2-mAb models than the bound protein: S1, chain H or chain L. The Zn-bound complexes also show lower RMSD values than their Zn-free counterparts at residues with at least one atom within 15 Å of Zn (Supplementary Information Figure S1). However, beyond this distance from Zn, the ACE2-S1 Zn-bound complex shows a greater RMSD than the Zn-free complex. To better understand the residues responsible for the changes in RMSD observed between the Zn-bound and Zn-free complexes, RMSF was computed for each residue within each protein (Fig. 1C–D). In the ACE2-S1 complex, RMSF of the Zn-bound complex is higher than the Zn-free complex at nearly every residue within ACE2 and S1 (Fig. 1C). In the ACE2-mAb complex, the Zn-bound complex exhibits lower RMSF at nearly every residue within ACE2. However, the difference in RMSF in the mAb chains H and L appears greatest in the first half of the residues, which are also more proximal to ACE2. The first half of residues within chains H and L of the mAb have greater RMSF in the Zn-free structure, while this difference is diminished or even reversed in the second half of each protein (Fig. 1D).

Fig. 1
figure 1

 A, B Average RMSD (left) of the interacting residues (IR) and the whole protein complex, both with and without zinc. Average RMSF (right) of each protein within the complex. These values were computed for both the ACE2-S1 complex (A), and the ACE2-mAb complex (B). Data represents mean ± s.e.m. C, D RMSF for each residue within the ACE2-S1 complex (C), and the ACE2-mAb complex (D). *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001 as determined by an unpaired t-test

Full size image

To investigate why the ACE2-S1 Zn-bound complex was more stable at the IR yet somewhat less stable overall, the time-based secondary structures of each complex were examined. A few differences are noticeable between the Zn-bound and Zn-free ACE2-S1 complexes (SI Figure S2). Unstable turns near ACE2 residue 107 in the Zn-bound complex are replaced by stable α-helices within the Zn-free complex. A similar change is also seen near ACE2 residue 448, where the Zn-bound complex contains unstable coils or 3–10 helices and the Zn-free complex contains stable α-helices. The secondary structures of the IR are of particular interest as it could potentially explain the deviations observed in Fig. 1. The secondary structures of IR within ACE2 show little variation between the complexes, although residue 355 of ACE2 is found to spend more time in stable β-sheet form in the Zn-bound complex (Fig. 2). Interestingly, in the ACE2-S1 system, more differences are noticed within the S1 protein’s secondary structure. In particular, S1 protein’s residues 449 and 496 appear to change from stable β-sheets in the Zn-bound complex to unstable turns and coils in the Zn-free complex (Fig. 2). In contrast, S1 residue 417 in the Zn-bound complex appears to fluctuate between stable α-helices and unstable turns and coils, while the residue remains as a stable α-helix in the Zn-free complex (Fig. 2). Some minor differences are observed in the secondary structure of the ACE2-mAb IR (Fig. 2). In the ACE2-mAb IR, residue 98 of chain H switches from turns in the Zn-bound complex to coils in the Zn-free complex. While this change in secondary structure is unlikely to directly alter RMSD, as both turns and coils are unstable, this shows that there are some variations in the ACE2-mAb interface.

Fig. 2
figure 2

Time-based changes in the secondary structures of the interacting residues (IR) in the ACE2-S1 complex (top) and ACE2-mAb complex (bottom), both in the Zn-bound (left) and Zn-free (right) complexes. In the ACE2-mAb complex, the inhibitor’s IR are split between chains H and L. The color at each horizontal position corresponds to a residue’s secondary structure at that time point in the 50 ns MD simulation

Full size image

Since we observed differences in the structural stability of the IR and total protein between Zn-bound and Zn-free complexes, we wanted to investigate if this is translated into a difference in binding free energy. Binding free energies were computed for the ACE2-S1 and ACE2-mAb complexes, both with and without Zn (Fig. 3). A range of internal dielectric constants were considered during binding energy calculation to account for charged and polar IR, which require higher internal dielectric constants for accurate computation. Since the ACE2-S1 interface contains a mixture of charged, polar, and nonpolar residues, the literature value for the binding free energy should be achievable with an internal dielectric constant less than 6. It was found that the binding free energies calculated with the internal dielectric constant of 2 align the most closely with the experimentally derived value of − 57.6 kcal/mol (Piplani et al. 2021). Since the mAb binds to the same site on ACE2, it is likely that the same internal dielectric constant can be used to compute ACE2-mAb binding free energy. For both the ACE2-S1 and ACE2-mAb systems, the Zn-bound complex has a stronger binding free energy than the Zn-free complex (Fig. 3). At an internal dielectric constant of 2, which aligns closely with the binding free energy literature values, these differences are − 3.26 and − 14.8 kcal/mol in the ACE2-S1 and ACE2-mAb systems respectively (Fig. 3).

Fig. 3
figure 3

Binding free energies of the ACE2-S1 complex (left) and the ACE2-mAb complex (right) computed using internal dielectric constants of 1, 2, 4, and 6. Both Zn-bound and Zn-free form of the complex were considered for computation. Data represents mean ± s.e.m. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001 as determined by an unpaired t-test

Full size image

Discussion

While Zn is expected to form strong connections to nearby residues within the protein and reduce movement of the overall complex, the results from MD simulation suggest that Zn may affect various protein-protein complexes differently. In comparison to the corresponding Zn-free complexes, the ACE2-S1 Zn-bound complex shows greater RMSD and lower stability while the ACE2-mAb Zn-bound complex shows an opposite behavior (Fig. 1A, B). In the ACE-S1 model, this difference in overall protein stability appears to be spread across the residues, as the RMSF of the Zn-bound complex appears greater at nearly all residues within ACE2 and S1 (Fig. 1C). Again, this finding differs from the ACE2-mAb complex, in which the greatest difference between the Zn-bound and Zn-free complexes is observed in the first half of the residues in chains H and L, which is more proximal to ACE2. Since this mAb was designed in-vivo against a Zn-containing ACE2, it is possible that the germinal center (GC) reaction has selected for an antigen binding domain that is extremely sensitive to the structural features of ACE2 (Du et al. 2021). Therefore, even small deviations in ACE2 structure from the exclusion of Zinc may have implications on the mAb’s binding configuration.

Although the exclusion of Zn can either increase or decrease the overall protein stability based on the given system, both Zn-bound complexes show greater stability in their IR compared to their Zn-free counterparts. In both the ACE2-S1 and ACE2-mAb models, the Zn-bound complex has IR with lower RMSD value than the corresponding Zn-free structure (Fig. 1A, B). In the ACE2-S1 model, some of these variations may be attributed to differences in the IR secondary structure. While the secondary structures for the ACE2 IR demonstrate only minor differences, the secondary structures of the S1 IR are found to be more stable in the Zn-bound complex (Fig. 2). For example, IR 449 and 496 in S1 change from more stable β-sheets in the Zn-bound complex to unstable turns and coils in the Zn-free complex. While Zn is reported to have distant effects (Sinha et al. 2016), at this time it is not fully understood why S1 residues located further away undergo changes even though much less changes are observed in the closer ACE2 residues. It is possible that the S1 residues are more susceptible to electrostatic and water-mediated hydrogen-bond interactions induced by the presence/absence of Zn, however, the reason for such susceptibility is yet to be understood.

Since differences in structural stability were observed, the next step was to check for differences in a common readout of inhibitor design: binding free energy. While computing the binding free energy, a variety of internal dielectric constants were used to account for charged and polar residues in the binding interfaces (Adekoya et al. 2006; Hou et al. 2011). Since S1 and mAb bind to the same region on ACE2, it is likely that the interface in both systems will contain a similar polarity and charge. Therefore, the optimal dielectric constant for the ACE2-S1 system can be applied to the ACE2-mAb system to estimate the binding free energy. Based on the literature value for the binding free energy of ACE2-S1, − 57.6 kcal/mol (Piplani et al. 2021) using an internal dielectric constant of two yields the most accurate results (Fig. 3). At this internal dielectric constant, there is a statistically significant difference in the binding free energy between the Zn-bound and Zn-free complexes in both the ACE2-S1 and ACE2-mAb systems. Interestingly the magnitude of the difference in binding free energy differs between the ACE2-S1 and ACE2-mAb systems, with differences of − 3.26 and − 14.8 kcal/mol respectively.

It is possible that the exclusion of Zn from the computational model affects the ACE2-mAb system to a greater degree than the ACE2-S1 system due to differences in the way that S1 and mAb were derived. In the GC reaction, there is a selective advantage for B cells which generate the B cell receptor (mAb) with the greatest antigen affinity. In contrast, while antigen specificity is important for SARS-CoV-2, there are many other selective pressures that this virus faces, which have even allowed newer strains to generate S1 domains with decreased affinity for ACE2 (Wu et al. 2022). Therefore, changes in protein structure, such as inclusion/exclusion of Zn, would likely have greatest impact on the binding of selectively designed proteins, such as mAb.

Conclusion

Designing inhibitor candidates for Zn-metalloproteins are quite challenging due to the metal’s complex coordination geometry and electron structure. Our results show that the exclusion of Zn from the ACE2-S1 and ACE2-mAb systems can influence structural stability and impacts the calculated value of binding free energy. Although coordination to Zn is expected to stabilize nearby residues, its inclusion can affect overall protein stability differently based on the system of interest. However, inclusion of Zn in both systems stabilized IR, which may explain observed improvements in binding free energy. The calculated binding free energies for both the ACE2-S1 and ACE2-mAb systems are stronger in the Zn-bound complexes than in the corresponding Zn-free complexes. It is possible that the difference in binding free energy is greatest in the ACE2-mAb system due to its generation from a GC reaction, which could make the mAb particularly sensitive to variations in the ACE2 structure. Since binding free energy is often used to score the efficacy of inhibitors, relatively slight differences in binding free energy values between the Zn-bound and Zn-free complex may lead to the inaccurate evaluation of the candidate inhibitors. This could also have implications for the design of inhibitors to other metalloproteins, which may also experience differences in the calculated binding free energy values from Zn exclusion.