Molecular Dynamics of the Neutralizing Potential of Dimeric IgA Antibodies against SARS-CoV-2 Antigen

Abstract

The emergence of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) stands for being the most serious epidemic (so far) of the 21^st century. However, only a few computational studies have investigated the molecular mechanisms underlying the neutralization of the spike protein by antibodies of different classes. Hence, bioinformatic methods were employed to unravel the factors contributing to the remarkable neutralization capacity exhibited by specific antibodies. Initially, crystallographic structures of IgA monomeric / dimeric, IgG, and IgM antibodies binding with the receptor-binding domain region of the SARS-CoV-2 spike protein were retrieved. Subsequently, rigid molecular docking and molecular dynamic simulations were performed over 100 ns with explicit water solvation. Lastly, an energy decomposition was conducted to estimate the binding affinity using the last frames from molecular dynamics. The results revealed a higher binding affinity for both monomeric and dimeric forms of IgA antibodies against the spike protein. Additionally, a greater number of hydrogen bonds were observed during their interaction with the spike protein, as well as greater structural instability along the time and especially a more thermodynamically favorable interaction affinity. In this way, the research contributes a small piece to the complex puzzle of understanding the humoral immune response induced by the SARS-CoV-2 virus.

Keywords: SARS-CoV-2; mucosal immunity; IgA antibodies; molecular dynamics

Introduction

The first cases of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection were reported in December 2019, in Wuhan, Hubei Province, China. This virus has caused a pandemic of a proportion not seen since the Spanish Flu pandemic of 1918, and, unfortunately, humanity was not prepared to deal with this new pandemic.^{1 2 3} In Brazil, as of October 24^th, 2023, more than 37.9 million cases have been confirmed and reported, with the live loss of around 706,000 lives.⁴ The SARS-CoV-2 virus encodes a total of four structural proteins, among which the most important in cell invasion are the spike glycoprotein (S) and the nucleocapsid protein (N).⁵, ⁶, ⁷

The (S) protein has a region called the receptor-binding domain (RBD), which interacts most strongly with the ACE2 receptor and is recognized by neutralizing antibodies. In addition, it is the main target in the development of vaccines, therapeutics, and molecular diagnostics..², ⁵, ⁸ Vaccination played a central and critical role in the fight against the Covid-19 (coronavirus disease 2019) pandemic by inducing an immunological memory response that saved thousands of lives, a fact had already been verified by using mathematical modelling.⁹, ¹⁰ Despite all the emergent variants of concern, such as P.1 / gamma that emerged in the city of Manaus (State of Amazonas, Brazil),¹¹ or the most recent Omicron variant, the critical importance of cellular immunity, where CD4+ and CD8+ T cells normally recognized the protein sequence spike, even with all it changes.¹, ¹², ¹³, ¹⁴

In the humoral response to SARS-CoV-2 infection, the literature has already established the importance of neutralizing antibodies, such as IgG, and this response is also closely correlated with the severity of the disease in some patients.¹⁵ Among the most important antibodies, there is the IgA class, which prevails in the human body initial response to the SARS-CoV-2 infection, compared to IgG and IgM concentrations.¹⁶, ¹⁷ Furthermore, IgA antibodies show greater virus neutralization than the IgG class. Thus, secretory IgA is predominant in the protection of the respiratory mucosa as a line of defense against pathogens.¹⁸ In addition, although IgG antibodies play a key role in virus neutralization, there is little information in the literature about the IgA-secreting antibodies that are predominant in the early stages of viral infection. Monomeric IgA antibodies for SARS-CoV-2 are found to be twice as potent as IgG, corroborated by experimental estudies.¹⁹, ²⁰, ²¹ One group of researchers analyzed the serological profile of patients with different levels of severity of Covid-19. In the results, it was noticed that the serum levels of IgA in patients with the severe form of the disease were more significant compared to patients with mild or moderate cases.², ²¹

In one study, the adenovirus vaccine commercialized by AstraZeneca (AZD1222) was administered to hamsters intranasally, and it was noticed that in the peripheral blood of the animals there were higher titers of neutralizing antibodies when compared to intramuscular vaccines.²⁰ In the enzyme-linked immunosorbent assay (ELISA), the IgA class of antibodies showed higher titers compared to IgG antibodies obtained from the peripheral blood of patients who died from Covid-19. Furthermore, there was a significant permanence of IgA antibodies for several days in patients who recovered from the acute infection.²² In a study carried out with indigenous populations in the city of Manaus, there were higher success rates for IgA compared to IgG when used as markers of SARS-CoV-2 infection. Furthermore, patients still had high levels of IgA even 4 weeks after first reporting their symptoms.¹⁶, ²², ²³

There have been some computational studies that sought to elucidate the SARS-CoV-2 infection at the molecular level, in particular, the role of the main protease (Mpro).⁵, ²³ However, theoretical simulations with different classes of antibodies neutralizing the spike protein are still scarce in the literature. In light of this, we used bioinformatics techniques such as molecular docking, molecular dynamics (MD) and the estimation of the free energy of interaction (∆G_binding) to examine the molecular antigen-antibody interaction. Finally, the Monte Carlo optimization method was used not only to find possible mutations that would increase the antibody-antigen interaction of SARS-CoV-2, but also to calculate theoretical emerging mutations in the future.

Methodology

Preparation of crystallographic structures

All the crystallographic structures utilized in this research were obtained from the RSCB Protein Data Bank (PDB).²⁴, ²⁵ However, the original files obtained from the PDB database were not totally suitable for computer simulations; thus, prior preparation was necessary. The structures were prepared using the “Protein Preparation Wizard” module integrated into the Schrödinger Maestro 2021-2 software.²⁶ During this process, hydrogen atoms that were not captured by the crystallography were added, partial charges were assigned and valence problems were corrected. In addition, to further refine the structures and eliminate steric clashes, structural minimization using the OPLS4 force field was performed.²⁷ The hydrogen atoms were optimized for a pH of 7.4, with the convergence criterion set to a root-mean square deviation (RMSD) lower than 0.30 Å. Due to computational constraints, the initial simulations focused only on the RBD of the spike protein. Among the simulated antibodies, only the IgG class was in complex with the spike protein. The specific antibody chosen for this work (PDB ID: 7BWJ) was found in the peripheral blood of patients affected by the severe form of Covid-19 with a resolution of 2.85 Å.²² In addition, it was also used another crystallography for the IgG antibody (PDB ID: 7BZ5) with a higher resolution of 1.84 Å, used at affinity maturation.²⁸ For the monomeric form of IgA1-Fc antibody, the crystallographic structure (PDB ID: 1OW0),¹⁸ obtained without being complexed with any antigen, was used. In addition, the IgM class antibody (PDB ID: 2AGJ)²⁹ and the dimeric IgA (PDB ID: 6UE7)³⁰ were chosen for this research. The available crystallographies did not specify the subclass of each antibody; therefore, the IgA, IgG, and IgM classes were analyzed in general.

Parameters of the molecular docking simulations

Due to the absence of available complexes of IgA and IgM antibodies with SARS-CoV-2 structural proteins, molecular docking simulations to explore how these antibodies bind to the spike protein was firstly performed. The “Protein-Protein Docking” module with the PIPER algorithm was utilized for all docking protocols.³¹ In the “Standard” mode of the software, the antibody was allowed 70,000 rotations, generating 30 possible conformation solutions. The conformation with the lowest potential energy was selected as the final docking pose. To quantify the affinity of the antibody-antigen complexes, molecular docking simulations were performed using three different docking tools: PatchDock,³² ClusPro,³³ and HDock.³⁴ Furthermore, to ensure the reliability of the interactions, molecular docking simulations were initiated from the last frame of molecular dynamics. The structural stability of all antibody-antigen complexes was then analyzed through molecular dynamic simulations.

Molecular dynamics system preparation

Initially, the complexes with best conformation of molecular docking results were retrieved. Then, molecular dynamic (MD) simulations were performed to analyze the most promising mutations that enhance antibody-antigen interaction based on Monte Carlo results. Molecular dynamic simulations with the neutralizing antibody IgG B38 (PDB ID: 7BZ5) isolated from the convalescent serum of a patient with COVID-19 were performed and,²⁸ the selectivity of antibodies was analyzed by also simulating IgA1-Fc (PDB ID: 1OW0).¹⁸ The Desmond algorithm³⁵ integrated into the Schrödinger Maestro 2021-2 software²⁶ was used for all the simulations. The systems were initially prepared using the “Preparation Wizard” module, whereby previously existing water molecules were removed, and hydrogen atoms were added and the protonation state assigned to pH ca. 7.4. Then, the “System Builder” function was used to build the TIP3P³⁶ cubic solvation box (Figure 1) for explicit water molecules and applying periodic boundary conditions while the system volume was bounded in a 10 Å box far from the complex (buffer distance). The exact dimensions were automatically determined by the minimize volume function. The systems were neutralized by adding Na⁺ and Cl⁻ ions at a concentration of 0.15 M. All systems were prepared and simulated using the OPLS4 force field. The simulated antibody-antigen systems were IgG B38 antibody (110,089 atoms with 33,474 water molecules), monomeric IgA1-Fc (71,134 atoms with 20,556 water molecules), dimeric IgA (351,810 atoms with 111,491 water molecules), IgM (107,936 atoms with 32,693 water molecules).

Figure 1
Solvation box created with “System Builder” function in Schrödinger Maestro 2021-2 software²⁶ antigen RBD with IgA-Fcl antibody (PDB ID: 1OW0) was inserted.

The system was minimized using the Desmond algorithm over 2 ns. Before beginning the simulations, the system was relaxed through equilibration. Finally, the simulation trajectory files were obtained in the isobaric-isothermal ensemble (NPT) ensemble at a temperature of 300.0 K and a pressure of 1 atm. A Nosé-Hoover thermostat³⁷ was employed to keep the temperature constant with a relaxation time of 1 ps and a Martyna-Tobias-Klein barostat³⁸ was used to maintain constant pressure with a relaxation time of 2 ps and isotropic coupling. The integration time in 2 fs was calculated using the RESPA algorithm, in which all the bonds surrounding the hydrogen atom were constrained using the M-SHAKE algorithm. Long-range electrostatic forces were calculated using the particle-mesh-Ewald (PME) technique³⁹ with a cut-off distance of 9.0 Å. The production of the trajectories was carried out at 300 K in the NPT ensemble over 100 ns. In addition, throughout all the simulations the OPLS4 force field was adopted. The simulations were performed with acceleration of the graphics processing unit (GPU) of model NVIDIA GTX 1050 2 GB with 640 CUDA cores in step of 2 fs.

Analysis of the molecular dynamics

Throughout all the analyses of the molecular dynamics, Schrödinger Maestro 2021-2²⁶ for RMSD, root-mean square fluctuations (RMSF), radius of gyration (Rgyr) and solvent accessible surface area (SASA) were used. Frame 1 was considered as a reference for the calculations where the entire protein complex was analyzed. In addition, Python 3 programming language⁴⁰ was used to create all the graphs used in this research. The average values of the molecular dynamics were calculated only for the last 20 ns of the simulation. In addition, statistical analyses were performed using Student’s f-test implemented in the Python 3 programming language using the Scipy library.⁴⁰ A 95% confidence interval for statistically significant differences was considered (p < 0.05). After the molecular dynamics at 100 ns, the Molecular-Mechanics General Born Surface Area (MM/GBSA) energy decomposition was performed from the last 5 ns, from which the Gibbs free energy (∆G_binding) of the antibody-antigen interaction was estimated using the OPLS4 force field and solvation was described using the VSGB model.⁴¹ Lastly, the chemical interaction networks were constructed using the iCn3D web platform⁴² from complex file resulting from MD simulations.

Affinity maturation for neutralizing antibodies

This computational method allows one to understand how different mutations in the spike (S) protein affect the binding affinity with antibodies, whether benefitting the virus by better understanding its evolutionary process, as well as the human being by increasing the potential of recognizing by neutralizing antibodies. In order to maximize the binding affinity and stability, the ACE2-RBD complex (PDB ID: 6M0J)³ and the antibody-antigen complexes for IgG B38 (PDB ID: 7BZ5)²⁸ and IgAl-Fc (PDB ID: 1OW0) were chosen.¹⁵ For this, the “Affinity Maturation” function integrated into the “Residue Scanning” module⁴³ within Schrödinger Maestro 2021-2²⁶ was used. In addition, MAESTROweb platform⁴⁴ was also employed for all analysis. This approach allowed us to conduct a Monte Carlo optimization search, which enabled the identification of specific mutations that enhance the binding affinity of the complexes. The property that was optimized was interaction affinity. This method was chosen because is computationally viable since it is not necessary to test the thousands of possible combinations of mutations. In order to make the simulations computationally viable, the solvent was implicitly considered using the MM-GBSA method in the Prime algorithm with the OPLS4 force field. This method was also fundamental for calculating all the energy values.

However, to remove possible steric conflicts, a structural minimization in the Prime algorithm was performed on the amino acids close to the mutation site within a cutoff of 5.0 Å. The total number of possible combinations for the 194 amino acids in the RBD was 4,074 mutations out of a total of 19 amino acids. At each step, one residue was randomly mutated. The initial seed was set to zero, while the temperature in all simulations was 300.0 K. The acceptance criteria for the next step occurred when the affinity / stability variation was greater than 30.0 kcal mol⁻¹. The process was repeated until the maximum number of steps was 200. Meanwhile, the maximum value for the number of simultaneous mutations was set to 16. These analyses were all performed with a physiological pH of 7.4.

Results and Discussion

Molecular dynamics of IgA and IgG antibodies

From the results of the molecular dynamics, the structural stability was quantified for the 100 ns of simulation (Figure 2 and Table 1). Firstly, as consequence of the un-specified crystallography subclass, the results more broadly for the IgG, IgM and IgA classes were analyzed and discussed. Thus, it was possible to realize the from the simulations that only RMSD (p < 0.05) and Rgyr (p < 0.05) results support the hypothesis of greater stability of monomeric IgA compared to IgG. However, the values of RMSF (p < 0.05) and SASA (p < 0.05) negate this hypothesis. A more stable system indicates greater formation of intermolecular bonds, which enhances antibody neutralization. Thus, the analysis of the chemical interactions formed in the antibody-antigen complex (Figure 3) indicates a greater neutralization of the antigen by the IgG type, as opposed to the monomeric form of IgA, because of a greater number of hydrogen bonds. Nevertheless, it nited that the dimeric form of the IgA antibody has greater neutralization potential, which is corroborated by the highest number of formed hydrogen bonds when, compared to the monomeric form, but also when compared to the IgG type.

Figure 2
(a) Comparison of RMSD values along 100 ns against the complex between RBD and antibodies IgA1 (monomer), IgG B38 and IgM; (b) comparison of RMSF fluctuation for the complex; (c) comparison of radius of gyration values for measurement of system compaction; and (d) comparison of SASA values for accessibility to solvent.

Table 1
Average values of RMSD, RMSF, Rgyr and SASA. Results for molecular dynamics simulations at 100 ns in Desmond software for the complex between RBD and two antibodies (IgG and IgA types)

Antibody	RMSD / Å	RMSF / Å	Rgyr / Å	SASA / nm2
IgG P2B-2F6	2.196 ± 0.328	1.281 ± 0.495	33.611 ± 0.197	292.05 ± 2.59
IgG B38	5.605 ± 1.027	1.690 ± 0.813	31.976 ± 0.19	283.18 ± 3.73
IgA1-Fc (monomer)	4.972 ± 0.921	2.105 ± 0.961	26.991 ± 0.316	301.97 ± 5.55
IgA (dimer)	31.697 ± 5.168	6.706 ± 5.456	43.890 ± 0.991	546.09 ± 4.28
IgM	7.152 ± 1.540	3.585 ± 1.314	33.086 ± 0.566	301.93 ± 4.17

• ^aRMSD (root mean square deviation) is used to measure the difference between the protein initial structure and each conformation over the course of the simulation; RMSF (root-mean square fluctuation) analysis is used to understand the individual contribution of each amino acid; Rgyr (radius of gyration) measures the degree of compaction of the system, and is also a measure of stability; SASA (solvent accessible surface area) is defined as the surface area of a protein that interacts with solvent molecules and which are most hydrophobic; summarizing, the stability of the protein over the simulation could be quantified by all these parameters.

Figure 3
(a) RMSD fluctuations analysis of the simulated IgA dimer antibody monomer over 100 ns in comparison with IgA monomer antibody; (b) RMSF analysis results; (c) radius of gyration (Rgyr) of system; (d) SASA values along the last 20 ns of simulation. All plots were generated with Python 3 / Matplotlib programming language.⁴⁰

The advantage of the RMSF analysis is that it is possible to know how each amino acid fluctuated individually in the antibody-antigen interaction. Thus, it can be noted that it is only after a certain region of amino acids that RMSF fluctuations become significant in the dimeric IgA (Figure 4). This observation is notable throughout the RMSF analyses, as it is possible to visualize the individual contribution of each amino acid. Thus, the RMSF chart was our option, as it seems to be the only one that reveals something fundamental in the interpretation of the IgA dimer. The rest of the graphs only show large peaks of fluctuations for the IgA antibody. Although the results of the molecular dynamics for the dimeric form of the IgA antibody seem inconsistent due to the very high structural instability over time, it is possible to notice that this behavior was repeated even after the simulations were redone. Consequently, it seems that high instability seems to be inherent to the phenomenon of spike neutralization by dimeric IgA antibodies, something that does not occur with monomeric IgA and other classes of antibodies.

Figure 4
Visualization of the last frame of the molecular dynamics for 100 ns of (a) monomeric IgA antibody (PDB ID: 1OW0) and (b) IgG antibody (PDB ID: 7BZ5) neutralizing the RBD of the spike protein. In addition, the RMSD result of the structural alignment between the last frame in relation to the first frame is presented; in red: spike protein, and in dark and light green heavy and light chains of the antibody, respectively.

Consequently, comparing only the first set (chain A and B) of amino acids from the dimeric IgA in complex with the spike RBD (chain E), it is possible to notice that the first fluctuations almost coincide with the monomeric IgA. However, the fluctuations start to become extremely high when analyzing precisely the second set of antibody chains and RBD. Therefore, the dimer induces substantial conformational changes in the second set of chains (F and G chain) that make up the IgA antibody. Therefore, it appears that the existence of the dimer is the main cause of greater fluctuations. The dimer consequently appears to lead to the appearance of a non-existent conformation in the monomer, perhaps an even more relaxed form of the IgA antibody.

The considerable number of fluctuations observed in the dimeric IgA highlights its enhanced flexibility and structural mobility during the interaction with RBD of the virus. Interestingly, the dimeric form of IgA demonstrated remarkable instability and flexibility when compared to the monomeric form, leading to a larger ∆S (entropy variation) and potentially rendering the value of ∆G (Gibbs free energy variation) more favorable for spike neutralization. The increased instability of the system suggests the inactivation of the spike protein, a key indicator of heightened neutralization by antibodies. These findings underscore the importance of antibody flexibility in facilitating effective viral neutralization and may contribute to the development of novel therapeutic strategies against SARS-CoV-2.

It can be observed that the average RMSD value of the molecular dynamics for the IgG antibody, RMSD ca. 2.20 ± 0.33 Å, was lower than that of the IgA antibody, RMSD ca. 4.97 ± 0.92 Å, with (p < 0.05), indicating greater structural stability in the antibody-antigen complex. This hypothesis is also reinforced by the RMSD peak value, which for the IgA antibody in its monomeric form was 6.75 Å, while, for the IgG class, the maximum value was 3.61 Å.

It is important to note that the index in the RMSF analyses (Figure 5) does not coincide with the position of the respective amino acid, since the crystallographic structures have missing residues. The RMSF value for the dimeric form of the IgA antibody increased sharply to 16.02 Å from the index position 646 equivalent to the amino acid Val524 of the spike protein, region from which the second IgA monomer also have very high values of atomic fluctuation. Meanwhile, the monomeric form of IgA peaked at a much lower RMSF, with a peak value of 7.34 Å at index position 418, which is equivalent to amino acid Thr333 of the spike protein.

Figure 5
RMSF plot as function of amino acid position for the dimeric form of the IgA antibody (PDB ID: 6UE7). The structure snapshots correspond to the last frame of the molecular dynamic simulation at 100 ns.

Free energy landscape

Firstly, Figure 6 shows the free energy landscape (FEL) for temperature (T) = 300 K, with smaller values corresponding to energetically favorable conformations. The panorama shows a clear separation between two wells for the case of the IgA dimer. One larger, with RMSD_min ca. 47.2 Å and another with RMSD_min ca. 44.9 Å. In the case of monomeric IgA, there is a well with RMSD_min ca. 5.5 Å and Rgyr greater than 44 Å and a region with Rgyr less than 44 Å with non-prominent local minima. For IgG, the single well has RMSD_min ca. 6.1 Å and Rgyr_min ca. 31.8 Å.

Figure 6
Free energy landscape as a function of RMSD and radius of gyration showing the energetically favored conformations. (a) IgA dimer energy surface; (b) IgA dimer main conformations; (c) IgA monomer energy surface; (d) IgA monomer main conformations; (e) IgG energy surface; and (f) IgG main conformation.

In Figures 6a and 6c, there is the surface map with the global configuration of energies. In Figures 6b, 6d and 6f, there is the three-dimensional map with the depth and spatial distribution of the wells. In Figure 6g, we have the three-dimensional structure corresponding to the energy wells for each case. In the case of Figures 6a and 6b, two wells are clearly observed, with the larger well having RMSD between 29.5-35.5 Å and an Rgyr value of less than 44 Å, with the smaller well having RMSD between 24-29.5 Å and an Rgyr value greater than 44 Å. At Figures 6c and 6d, a well with RMSD between 5-6 Å and Rgyr value greater than 27 Å with local minima for Rgyr intervals of lower than 27 Å is observed. In Figure 6f, we have a single well with a RMSD_min ca. 6.1 Å and Rgyr_mn ca. 31.8 Å.

A wider distribution along the RMSD axis can be seen for the dimeric IgA compared to the monomeric IgA and IgG. Thus, it is observed that the monomeric structure of IgA, as well as of IgG, did not change significantly when compared to dimeric IgA; moreover, both monomeric IgA and IgG have a greater compaction (smaller range of Rgyr values). The dimeric structure, as it has two well-defined valleys, demonstrates that the minimum free energy becomes unstable, and therefore has two predominant conformations. The monomeric structure is less unstable since the secondary minima are not located in deep valleys. The IgG structure, on the other hand, has a single well-defined valley which makes a predominant and stable conformation.

MM-GBSA energy decomposition and molecular docking

By decomposing the energy with the latest molecular dynamics structure frame, it was noted that the monomeric form of the IgA antibody actually provided more favorable values in the Gibbs free energy, with −207.293 kcal mol⁻¹. Secondly, we have the dimeric form of IgA (−199.120 kcal mol⁻¹), which showed better results to the IgG type (−133.805 kcal mol⁻¹). Finally, we had the IgM antibody with extremely low neutralizing potential, with −22.058 kcal mol⁻¹, although has a significative importance for first line immunological defense. The theoretical results are, therefore, in line with what has been reported experimentally, where dimeric IgA stands out with a much higher neutralization potential than the others (Table 2).⁹ We should note that these free energy values are relative, so we can only perform a qualitative comparison between antibodies consistent with the experimental results.

Table 2
Results of MM-GBSA energy (∆G_binding) decomposition method for antibody-antigen complexes. All the results correspond to the average of last 5 ns of the molecular dynamic simulations over 100 ns

Antibody	AGbinding / (kcal mol−1)
IgG B38	−133.805
IgA-Fc1 (monomer)	−207.293
IgA (dimer)	−199.120
IgM	−22.058

Despite the huge importance of molecular dynamic simulations, interpretations and analysis based only on this technique are very abstract, so we performed molecular docking to get a better idea of the change in affinity between the different antibodies (Table 2). From the molecular docking analysis (Table 3), the dimeric form of the SIgA antibody was absolutely superior in all the four tools adopted. In this structural form of the IgA antibody, as a consequence, we have a crucial response in the face of infection by the SARS-CoV-2 virus. In comparison, the IgG antibody showed the second highest affinity in three tools, Prodigy,³⁵ HDock³⁴ and PatchDock.³² Meanwhile the IgA type (monomer) was higher only in the ClusPro tool.⁴⁵ This is an indication that the IgG type has a greater neutralizing potential, but the result is completely inverted when compared with dimeric IgA.

Table 3
Molecular docking scores for the last frame in the molecular dynamic simulations of antigen-antibody complex at 100 ns

Antibody	Prodigy	HDock	PatchDock	ClusPro
IgG B38	−14.6	−386.00	12,616	−861.4
IgA1-Fc (monomer)	−10.7	−313.35	11,218	−928.7
IgA (dimer)	−17.8	−587.44	13,598	−1,458.9
IgM	−8.4	−244.94	12,512	−768.2

Diagram of intermolecular interactions

To gain insights into the neutralization potential of different antibody types against the spike protein, the chemical interactions formed in the antibody-antigen complex were quantified (Figure 7). First, the analysis was performed with Schrödinger Maestro 2021-2,²⁶ and revealed distinct patterns of hydrogen bonding in the complex for each antibody type. In the case of IgA (monomer), a total of 9 (nine) hydrogen bonds formed with the spike protein chain was observed, as quantified from the last frame of the molecular dynamic simulation. On the other hand, the IgG type exhibited a higher number of interactions, with a total of 14 (fourteen) hydrogen bonds identified in the complex. Notably, the dimeric form of the IgA antibody displayed a significant increase in interactions, with a total of 20 (twenty) hydrogen bonds established between the spike protein chain and the antibody chains. These findings shed light on the distinct binding modes and stability of the antibody-antigen complexes, and suggest that the dimeric form of IgA may possess a superior neutralizing capacity against the spike protein when compared to the monomeric IgA and IgG antibodies.

Figure 7
Diagram of chemical interactions formed in the last frame of the molecular dynamic simulations at 100 ns for the complex between the spike protein and the monomeric IgA antibody (PDB ID: 1OW0), in addition to the IgM antibody (PDB ID: 2AGJ). The purple arrows represent the hydrogen bonds formed. The chain E represents the spike protein, while the chains A, B, H and L represent the neutralizing antibody.

By comparing the interaction networks (Figure 8) formed between monomeric IgA and IgG with an approach of another software, ICn3d,⁴² it can better be visualized their molecular mechanisms of neutralization, corroborating or not the previously analysis. Thus, it was noted that, for IgA, 3 (three) salt bridges (in light blue) were formed, in addition to 4 (four) hydrogen bonds (in green). On the other hand, in view of the neutralization by the IgG class antibody, only 2 (two) hydrogen bonds and 2 (two) salt bridges were formed, in addition to the appearance of one π-stacking bond (in dark blue). Thus, it is noted that the IgA class showed greater intermolecular binding when neutralizing the RBD of the spike protein in comparison to the IgG antibody. The difference between the number of hydrogen bonds formed previously in comparison with ICn3d⁴² is due to the algorithm used in each software. However, the main purpose is just to confirm that IgA antibodies actually induce the formation of a greater number of hydrogen bonds. These interactions are also reflected in a more spontaneous value of free energy, as verified by the molecular dynamics. Finally, although interactions tend to change throughout molecular dynamics, we chose only the last frame for more consistent results.

Figure 8
Network of chemical interactions formed in the neutralization of the RBD of the spike protein by monomeric IgA antibody (PDB ID: 1OW0) and also by IgG antibody (PDB ID: 7BWJ).

Affinity maturation results

Among the analyses performed, the affinity maturation searches for maximized affinity that benefits human beings (Table S6 in the Supplementary Information (SI) section), which would increase the antibody-antigen interaction. Ultimately, the minimized affinity would benefit the virus as the interaction between antibody and antigen would become minimal. The results for minimized affinity on the antibody-antigen interaction reflect the worst possibility for the virus to be able to fully escape from the neutralizing antibodies. Nevertheless, predicting the viral behavior and its mutations is unfeasible with these techniques, and, therefore, it was only aimed to maximize the antibody-antigen affinity.

By using the Monte Carlo method, it was found a total of sixteen mutations in the antigen (Table S7 in the SI section) that resulted in a maximized affinity of −61.656 kcal mol⁻¹ between the IgG antibody B38 and the spike protein (PDB ID: 7BZ5). Among these, it was identified G339Q, T345R, F347C, S349C, I358L, N360D, Y380R, G416R, S438F, L455R, F490L, Y508I, L518M, A520V and P527V. On the other hand, in relation to the IgA1-Fc antibody (PDB ID: 1OW0), it was found that the maximized affinity with the antigen would be -82.343 kcal mol⁻¹. Thus, these mutations were: T333V, N370M, S373E, S399R, I402Q, D405W, Y421S, N422H, Y449L, S469N, E471T, P491Y, S494P, Y495R, F497H and F515R.

We realized that IgG, IgM and IgA antibodies showed different combinations of mutations that maximize the interaction affinity (Table S2 in the SI section). This is an indication of selectivity, in other words, a certain set of mutations will only benefit the neutralization of one antibody over another. Nevertheless, performing simulations for different classes of antibodies could contribute to obtaining a more selective response depending on the chosen mutations for the antigen.

Using the MAESTROweb platform,⁴⁴ it was possible to have a better idea about theoretical mutations that could emerge in the future, where was searched those caused a decrease in virus affinity for neutralizing antibodies as a function of the pH condition (Table S2 in the SI section). The choice of pH values was to simulate the gastric environment in the acidic condition at pH ca. 5.5, enterocytes at pH ca. 6.8, the pulmonary epithelium in the physiological / slightly alkaline condition at pH ca. 7.4, and, finally, the intestine in the basic condition with pH ca. 8.5. Interestingly, when analyzing the dependence on pH for the appearance of new mutations, the slightly acidic condition with a pH of 6.8 favored the appearance of less unstable mutations compared to the other conditions, for which we obtained ∆G ca. +5.383 kcal mol⁻¹, while, in alkaline condition, we obtained ∆G ca. +3.649 kcal mol⁻¹ and, in the physiological condition, the value of ∆G ca. +4.985 kcal mol⁻¹. This may explain why the acidic condition is more favorable for viral replication.

Limitations of the study

In highlight, the time interval of 100 ns was adopted for molecular dynamic simulation and the lack of replicates, although there is a convergence of fluctuations after the system reaches equilibrium. Although the application of the Monte Carlo method for SARS-CoV-2 is certainly the most promising result of this research, all the calculations of ∆G_binding were based on a static protein structure, instead of dynamic behavior. Due to the enormous advances in artificial intelligence (AI), especially the deep-learning approach, the use of tools such as AlphaFold for docking between proteins could better estimate how antibodies and antigens interact when there is no crystallographic structure of the complex, such as IgA antibodies of this research. Finally, further studies on the dynamics and avidity of mucosal IgA neutralizing antibodies are still necessary in order to better understand its mechanism in vivo.

Conclusions

Given the significance of mucosal immunity mediated by the IgA class, molecular dynamic simulations were performed to investigate its enhanced affinity in neutralizing the spike protein. The molecular dynamics results were quantified by RMSD, RMSF, Rgyr, and SASA values, as also binding free energy. In this way, leading to the conclusion regarding the primary mechanism of inhibition by the dimeric form of the IgA antibody class, that underscored significant conformational changes throughout the simulation. Thus, it was observed that IgA exhibited greater structural instability and formed a higher number of hydrogen bonds compared to other antibody classes. Furthermore, we proposed protein engineering with specific antigen mutations that could enhance susceptibility to neutralization by antibodies by increasing the binding affinity of interaction, which could be valuable in potential treatments if there is a new pandemic. In summary, this study suggests that IgA antibodies may have greater immunogenicity if futures vaccines against COVID-19 or other respiratory infections were administered via the nasal route where IgA class predominates.

Supplementary Information

All trajectories, topology and analysis files for molecular dynamics at 100 ns are available on the Zenodo platform (https://zenodo.org/records/7583111). Supplementary data are available free of charge at https://jbcs.sbq.org.br as PDF file.

Acknowledgments

We express our gratitude to all the researchers and other people who came together to help each other in the fight against the SARS-CoV-2 virus. We would like to thank National Council for Scientific and Technological Development (CNPq) for supplying this research with grant 128869/2021-6. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior, (CAPES, Brasil), Finance Code 001.

References

• ¹ Li, X.; Giorgi, E. E.; Marichannegowda, M. H.; Foley, B.; Xiao, C.; Kong, X.-P.; Chen, Y.; Gnanakaran, S.; Korber, B.; Gao, F.; Sci. Adv. 2020, 6, eabb9153. [Crossref]
  » Crossref
• ² Harrison, A. G.; Lin, T.; Wang, P.; Trends Immunol. 2020, 41, 1100. [Crossref]
  » Crossref
• ³ Lan, J.; Ge, J.; Yu, J.; Shan, S.; Zhou, H.; Fan, S.; Zhang, Q.; Shi, X.; Wang, Q.; Zhang, L.; Wang, X.; Nature 2020, 581, 215. [Crossref]
  » Crossref
• ⁴ Painel Coronavírus SUS, https://covid.saude.gov.br/, accessed in October 10, 2023.
  » https://covid.saude.gov.br/
• ⁵ Arya, R.; Kumari, S.; Pandey, B.; Mistry, H.; Bihani, S. C.; Das, A.; Prashar, V.; Gupta, G. D.; Panicker, L.; Kumar, M.; J. Mol. Biol. 2021, 433, 166725. [Crossref]
  » Crossref
• ⁶ Malone, B.; Urakova, N.; Snijder, E. J.; Campbell, E. A.; Nat. Rev. Mol. Cell Biol. 2022, 23, 21. [Crossref]
  » Crossref
• ⁷ V’kovski, P.; Kratzel, A.; Steiner, S.; Stalder, H.; Thiel, V.; Nat. Rev. Microbiol. 2021, 19, 155. [Crossref]
  » Crossref
• ⁸ Wrapp, D.; Wang, N.; Corbett, K. S.; Goldsmith, J. A.; Hsieh, C.-L.; Abiona, O.; Graham, B. S.; McLellan, J. S.; Science 2020, 367, 1260. [Crossref]
  » Crossref
• ⁹ Watson, O. J.; Barnsley, G.; Toor, J.; Hogan, A. B.; Winskill, P.; Ghani, A. C.; Lancet Infect. Dis. 2022, 22, 1293. [Crossref]
  » Crossref
• ¹⁰ Wang, M.-Y.; Zhao, R.; Gao, L.-J.; Gao, X.-F.; Wang, D.-P.; Cao, J.-M.; Front. Cell. Infect. Microbiol. 2020, 10, 587269. [Crossref]
  » Crossref
• ¹¹ Naveca, F. G.; Nascimento, V.; de Souza, V. C.; Corado, A. L.; Nascimento, F.; Silva, G.; Costa, Á.; Duarte, D.; Pessoa, K.; Mejía, M.; Brandão, M. J.; Jesus, M.; Gonçalves, L.; da Costa, C. F.; Sampaio, V.; Barros, D.; Silva, M.; Mattos, T.; Pontes, G.; Abdalla, L.; Santos, J. H.; Arantes, I.; Dezordi, F. Z.; Siqueira, M. M.; Wallau, G. L.; Resende, P. C.; Delatorre, E.; Gräf, T.; Bello, G.; Nat. Med. 2021, 27, 1230. [Crossref]
  » Crossref
• ¹² Tarke, A.; Sidney, J.; Methot, N.; Yu, E. D.; Zhang, Y.; Dan, J. M.; Goodwin, B.; Rubiro, P.; Sutherland, A.; Wang, E.; Frazier, A.; Ramirez, S. I.; Rawlings, S. A.; Smith, D. M.; da Silva Antunes, R.; Peters, B.; Scheuermann, R. H.; Weiskopf, D.; Crotty, S.; Grifoni, A.; Sette, A.; Cell Rep. Med. 2021, 2, 100355. [Crossref]
  » Crossref
• ¹³ Suzuki, R.; Yamasoba, D.; Kimura, I.; Wang, L.; Kishimoto, M.; Ito, J.; Morioka, Y.; Nao, N.; Nasser, H.; Uriu, K.; Kosugi, Y.; Tsuda, M.; Orba, Y.; Sasaki, M.; Shimizu, R.; Kawabata, R.; Yoshimatsu, K.; Asakura, H.; Nagashima, M.; Sadamasu, K.; Yoshimura, K.; Suganami, M.; Oide, A.; Chiba, M.; Ito, H.; Tamura, T.; Tsushima, K.; Kubo, H.; Ferdous, Z.; Mouri, H.; Iida, M.; Kasahara, K.; Tabata, K.; Ishizuka, M.; Shigeno, A.; Tokunaga, K.; Ozono, S.; Yoshida, I.; Nakagawa, S.; Wu, J.; Takahashi, M.; Kaneda, A.; Seki, M.; Fujiki, R.; Nawai, B. R.; Suzuki, Y.; Kashima, Y.; Abe, K.; Imamura, K.; Shirakawa, K.; Takaori-Kondo, A.; Kazuma, Y.; Nomura, R.; Horisawa, Y.; Nagata, K.; Kawai, Y.; Yanagida, Y.; Tashiro, Y.; Takahashi, O.; Kitazato, K.; Hasebe, H.; Motozono, C.; Toyoda, M.; Tan, T. S.; Ngare, I.; Ueno, T.; Saito, A.; Butlertanaka, E. P.; Tanaka, Y. L.; Morizako, N.; Sawa, H.; Ikeda, T.; Irie, T.; Matsuno, K.; Tanaka, S.; Fukuhara, T.; Sato, K.; Nature 2022, 603, 700. [Crossref]
  » Crossref
• ¹⁴ Wang, Q.; Guo, Y.; Iketani, S.; Nair, M. S.; Li, Z.; Mohri, H.; Wang, M.; Yu, J.; Bowen, A. D.; Chang, J. Y.; Shah, J. G.; Nguyen, N.; Chen, Z.; Meyers, K.; Yin, M. T.; Sobieszczyk, M. E.; Sheng, Z.; Huang, Y.; Liu, L.; Ho, D. D.; Nature 2022, 608, 603. [Crossref]
  » Crossref
• ¹⁵ Qi, H.; Liu, B.; Wang, X.; Zhang, L.; Nat. Immunol. 2022, 23, 1008. [Crossref]
  » Crossref
• ¹⁶ Sterlin, D.; Mathian, A.; Miyara, M.; Mohr, A.; Anna, F.; Claër, L.; Quentric, P.; Fadlallah, J.; Devilliers, H.; Ghillani, P.; Gunn, C.; Hockett, R.; Mudumba, S.; Guihot, A.; Luyt, C.-E.; Mayaux, J.; Beurton, A.; Fourati, S.; Bruel, T.; Schwartz, O.; Lacorte, J.-M.; Yssel, H.; Parizot, C.; Dorgham, K.; Charneau, P.; Amoura, Z.; Gorochov, G.; Sci. Transl. Med. 2021, 13, eabd2223. [Crossref]
  » Crossref
• ¹⁷ Wang, Z.; Lorenzi, J. C. C.; Muecksch, F.; Finkin, S.; Viant, C.; Gaebler, C.; Cipolla, M.; Hoffmann, H.-H.; Oliveira, T. Y.; Oren, D. A.; Ramos, V.; Nogueira, L.; Michailidis, E.; Robbiani, D. F.; Gazumyan, A.; Rice, C. M.; Hatziioannou, T.; Bieniasz, P. D.; Caskey, M.; Nussenzweig, M. C.; Sci. Transl. Med. 2021, 13, eabf1555. [Crossref]
  » Crossref
• ¹⁸ Herr, A. B.; Ballister, E. R.; Bjorkman, P. J.; Nature 2003, 423, 614. [Crossref]
  » Crossref
• ¹⁹ Ma, H.; Zeng, W.; He, H.; Zhao, D.; Jiang, D.; Zhou, P.; Cheng, L.; Li, Y.; Ma, X.; Jin, T.; Cell. Mol. Immunol. 2020, 17, 773. [Crossref]
  » Crossref
• ²⁰ Fischer, R. J.; van Doremalen, N.; Adney, D. R.; Yinda, C. K.; Port, J. R.; Holbrook, M. G.; Schulz, J. E.; Williamson, B. N.; Thomas, T.; Barbian, K.; Anzick, S. L.; Ricklefs, S.; Smith, B. J.; Long, D.; Martens, C.; Saturday, G.; de Wit, E.; Gilbert, S. C.; Lambe, T.; Munster, V. J.; Nat. Commun. 2021, 12, 5868. [Crossref]
  » Crossref
• ²¹ de-Oliveira-Pinto, L. M.; Fiestas Solorzano, V. E.; de Lourdes Martins, M.; Fernandes-Santos, C.; Damasco, P. H.; de Siqueira, M. A. M. T.; Dias, H. G.; Pauvolid-Corrêa, A.; Damasco, P. V.; de Azeredo, E. L.; Viruses 2022, 14, 455. [Crossref]
  » Crossref
• ²² Ju, B.; Zhang, Q.; Ge, J.; Wang, R.; Sun, J.; Ge, X.; Yu, J.; Shan, S.; Zhou, B.; Song, S.; Tang, X.; Yu, J.; Lan, J.; Yuan, J.; Wang, H.; Zhao, J.; Zhang, S.; Wang, Y.; Shi, X.; Liu, L.; Zhao, J.; Wang, X.; Zhang, Z.; Zhang, L.; Nature 2020, 584, 115. [Crossref]
  » Crossref
• ²³ Biembengut, Í. V.; de Souza, T. A. C. B.; Mem. Inst. Oswaldo Cruz 2020, 115, e190401. [Crossref]
  » Crossref
• ²⁴ RSCB Protein Data Bank (PDB), https://www.rcsb.org/,accessed in July 2024.
  » https://www.rcsb.org/
• ²⁵ Burley, S. K.; Bhikadiya, C.; Bi, C.; Bittrich, S.; Chen, L.; Crichlow, G. V; Christie, C. H.; Dalenberg, K.; Di Costanzo, L.; Duarte, J. M.; Dutta, S.; Feng, Z.; Ganesan, S.; Goodsell, D. S.; Ghosh, S.; Green, R. K.; Guranovic, V.; Guzenko, D.; Hudson, B. P.; Lawson, C. L.; Liang, Y.; Lowe, R.; Namkoong, H.; Peisach, E.; Persikova, I.; Randle, C.; Rose, A.; Rose, Y.; Sali, A.; Segura, J.; Sekharan, M.; Shao, C.; Tao, Y.-P.; Voigt, M.; Westbrook, J. D.; Young, J. Y.; Zardecki, C.; Zhuravleva, M.; Nucleic. Acids Res. 2021, 49, D437. [Crossref]
  » Crossref
• ²⁶ Maestro-Desmond Interoperability Tools; SchrödingerRelease 2021-2: Desmond Molecular Dynamics System; Schrödinger Inc., New York, NY, 2021.
• ²⁷ Lu, C.; Wu, C.; Ghoreishi, D.; Chen, W.; Wang, L.; Damm, W.; Ross, G. A.; Dahlgren, M. K.; Russell, E.; Von Bargen, C. D.; Abel, R.; Friesner, R. A.; Harder, E. D.; J. Chem. Theory Comput. 2021, 17, 4291. [Crossref]
  » Crossref
• ²⁸ Wu, Y.; Wang, F.; Shen, C.; Peng, W.; Li, D.; Zhao, C.; Li, Z.; Li, S.; Bi, Y.; Yang, Y.; Gong, Y.; Xiao, H.; Fan, Z.; Tan, S.; Wu, G.; Tan, W.; Lu, X.; Fan, C.; Wang, Q.; Liu, Y.; Zhang, C.; Qi, J.; Gao, G. F.; Gao, F.; Liu, L.; Science 2020, 568, 1274. [Crossref]
  » Crossref
• ²⁹ Ramsland, P. A.; Terzyan, S. S.; Cloud, G.; Bourne, C. R.; Farrugia, W.; Tribbick, G.; Geysen, H. M.; Moomaw, C. R.; Slaughter, C. A.; Edmundson, A. B.; Biochem. J. 2006, 595, 473. [Crossref]
  » Crossref
• ³⁰ Kumar, N.; Arthur, C. P.; Ciferri, C.; Matsumoto, M. L; Science 2020, 567, 1008. [Crossref]
  » Crossref
• ³¹ Kozakov, D.; Brenke, R.; Comeau, S. R.; Vajda, S.; Proteins Struct. Funct. Bioinf. 2006, 65, 392. [Crossref]
  » Crossref
• ³² Schneidman-Duhovny, D.; Inbar, Y.; Nussinov, R.; Wolfson, H. J.; Nucleic Acids Res. 2005, 55, W363. [Crossref]
  » Crossref
• ³³ Kozakov, D.; Hall, D. R.; Xia, B.; Porter, K. A.; Padhorny, D.; Yueh, C.; Beglov, D.; Vajda, S.; Nat. Protoc. 2017, 12, 255. [Crossref]
  » Crossref
• ³⁴ Yan, Y.; Zhang, D.; Zhou, P.; Li, B.; Huang, S.-Y.; Nucleic Acids Res. 2017, 45, W365. [Crossref]
  » Crossref
• ³⁵ Shaw, D. E.; Maragakis, P.; Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Eastwood, M. P.; Bank, J. A.; Jumper, J. M.; Salmon, J. K.; Shan, Y.; Wriggers, W.; Science 2010, 550, 341. [Crossref]
  » Crossref
• ³⁶ Mark, P.; Nilsson, L.; J. Phys. Chem. A 2001, 105, 9954. [Crossref]
  » Crossref
• ³⁷ Evans, D. J.; Holian, B. L.; J. Chem. Phys. 1985, 85, 4069. [Crossref]
  » Crossref
• ³⁸ Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L.; Mol. Phys. 1996, 87, 1117. [Crossref]
  » Crossref
• ³⁹ Petersen, H. G.; J. Chem. Phys. 1995, 105, 3668. [Crossref]
  » Crossref
• ⁴⁰ Van Rossum, G.; Drake Jr., F. L.; Python 5 Programming Language with Libraries Matplotlib and Scipy; Python Software Foundation, 2022.
• ⁴¹ Li, J.; Abel, R.; Zhu, K.; Cao, Y.; Zhao, S.; Friesner, R. A.; Proteins Struct. Funct. Bioinf. 2011, 79, 2794. [Crossref]
  » Crossref
• ⁴² Wang, J.; Youkharibache, P.; Zhang, D.; Lanczycki, C. J.; Geer, R. C.; Madej, T.; Phan, L.; Ward M.; Lu, S.; Marchler, G. H.; Wang, Y.; Bryant, S. H.; Lewis Y Geer, L. Y.; Marchler-Bauer, A.; Bioinformatics 2020, 56, 131 [Crossref] [PubMed]; Wang, J.; Youkharibache, P.; Marchler-Bauer, A.; Lanczycki, C.; Zhang, D.; Lu S.; Madej, T.; Marchler, G. H.; Cheng, T.; Chong, L. C.; Zhao, S.; Yang, K.;, Lin, J.; Cheng, Z.; Dunn, R.; Malkaram, S. A.; Tai, C.-H.; Enoma, D.; Busby, B.; Johnson, N. L.; Tabaro, F.; Song, G.; Ge, Y.; Front. Mol. Biosci. 2022, 9, 831740. [Crossref] [PubMed]
  » Crossref» PubMed» Crossref» PubMed
• ⁴³ Beard, H.; Cholleti, A.; Pearlman, D.; Sherman, W.; Loving, K. A.; PLoS One 2013, 8, e82849. [Crossref]
  » Crossref
• ⁴⁴ Hiebl-Flach J.; Laimer, J.; Laimer, J.; MAESTROweb; University of Salzburg, 2015, [Link] accessed in July 2024; Laimer, J.; Hiebl-Flach, J.; Lengauer, D.; Lackner, P.; Bioinformatics 2016, 52, 1414. [https://doi.org/10.1093/bioinformatics/btv769">Crossref]
  » Link
• ⁴⁵ Pontes, G. S.; de Melo Silva, J.; Pinheiro-Silva, R.; Barbosa, A. N.; Santos, L. C.; de Pádua Quirino Ramalho, A.; de Castro Alves, C. E.; da Silva, D. F.; de Oliveira, L. C.; da Costa, A. G.; Bruno, A. C.; Sci. Rep. 2021, 11, 17534. [Crossref]
  » Crossref

Edited by

• Editor handled this article: Paula Homem de Mello (Associate)

Publication Dates

• Publication in this collection
  23 Sept 2024
• Date of issue
  2025

History

• Received
  15 Feb 2024
• Accepted
  09 Aug 2024