In silico antibody engineering for SARS-CoV-2 detection

Graphical abstract


Introduction
Immunoglobulin G (IgG) antibodies, which play a key role regulating the human immune system [1], are amongst the most exquisitely designed and engineered molecules in Nature. Because of their exceptional bio-recognition elements, which exhibit high specificity and affinity for their cognate antigen, IgG antibodies are widely used in serological sensor devices (immunosensors) for detection of pathogens and toxins [2][3][4][5][6][7][8]. Within this context, recombinant technology, which easily allows antibodies genetic manipulation, is a valuable and robust tool for the fabrication of immunosensors.
IgG antibodies present a monomeric ''H2L2" structure, consisting of two identical heavy chains (H) and two identical light (L) chains. The molecular weight of H chains ($50 kDa) is approximately twice that of L chains ($25 kDa), the former ones being linked among them and to a L chain each by disulfide bonds. The quaternary structure is characterized by two identical halves that joint forming a Y-like shape (Scheme 1). Each H chain contains a variable domain (V H ) and three constant domains (CH 1 , CH 2 , CH 3 ), with an additional ''hinge region" between CH 1 and CH 2 . Besides, L chains involve variable domain (V L ) and a constant domain (C L ). The light chain associates with the V H and CH 1 domains to form the fragment antigen binding (Fab) arm (''Fab" = fragment antigen binding). The lower hinge region and the CH 2 /CH 3 domains form the fragment crystalline (Fc). The fragment composed of Fab regions, joined by the hinge region, is known as F (ab') 2 .
The performance of manufactured antibody-based sensors depends, among others, on the procedure used to immobilize antibody while maintaining its natural activity. Thus, although antibodies can be attached to the solid support in many different orientations, immobilization through the Fc region is crucial to leave the Fab exposed for recognition of the antigen, optimizing https://doi.org/10.1016/j.csbj.2021. 10 the sensitivity and the limit of detection of the sensor [9]. Within this context, different strategies have been reported to avoid random orientations, promoting the immobilization of the Fcspecific orientation (e.g. photon-assisted methods based on ultrashort UV laser pulses [10], controlled electric fields [11] and surface functionalization [9,12]).
Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) has triggered a global health crisis with high social impact through the COVID-19 disease [13]. The spike glycoprotein protein, which consists on a trimer with three monomers exhibiting identical primary structure, plays a key role in viral infection and pathogenesis [14]. SARS-CoV-2 infection undergoes a series of processes. The binding of the receptor-binding domain (RBD) to its receptor, angiotensin converting enzyme 2 (ACE2), to form an RBD/ACE2 complex is the first [15,16]. It triggers conformational changes in the spike protein, leading to membrane fusion mediated via others part of the spike [17,18]. This process culminates in viral entry into target cells. Within this context, the development of rapid and efficient immunosensors for early detection in the diagnosis of SARS-CoV-2 is highly desirable.
The design of neutralizing IgG-like antibodies that binds the RBD of the SARS-CoV-2 protein has been found to be a promising approach against COVID-19 [19][20][21][22]. Indeed, a huge amount of SARS-CoV-2 antibodies, obtained from COVID-19 patients, have been disclosed since the beginning of the pandemic [23]. In this work we focus on two of them, CR3022 and S309. CR3022 does not overlap and not compete with the ACE2 binding site when binds SARS-CoV-2 [24]. Although there is some controversy on the in vivo neutralizing capacity of CR3022 against SARS-CoV-2, the antibody is able to bind to the RBD, conferring in vivo protection [25][26][27]. The binding of the CR3022 to the epitope on RBD of SARS-CoV-2 spike was recently studied using structural Molecular Dynamics (MD) simulations [28]. On the other hand, S309 monoclonal antibody was recently identified as a therapeutic agent that potently neutralizes SARS-CoV-2 by recognizing an epitope that contains a fucosylated glycan at position N343 [29]. S309 binds the RBD of SARS-CoV-2 spike in open conformation without competing with ACE2 binding. Currently, an engineered antibody (VIR-7831) based on S309 is under clinical trial assessment in humans (i.e. NCT04545060) [30,31].
In this work, two IgG-like antibodies have been specifically engineered to detect SARS-CoV-2. After that, in silico studies have been conducted to examine their stability when immobilized on a solid gold surface, which is a usual substrate for virus immunosensors [32][33][34][35]. For this purpose, the crystal structure of the neutralizing antibody CR3022 Fab tightly bounded to SARS-CoV-2 RBD (Code in the Protein Data Bank, PDB: 6YLA) [25] and the crystal structure of the S309 neutralizing antibody Fab fragment in complex with SARS-CoV-2 RBD (Code in the PDB: 6WPT) [29] were combined with the crystal structure of IgG1 B12 antibody (Code in PDB: 1HZH) [36], which was used to inactivate the human immunodeficiency virus-1 (HIV-1). More specifically, the Fab fragment of CR3022 and the Fv fragment (i.e. the fragment formed by V H and V L , which are held together by non-covalent interactions; see Scheme 1) of S309 were used to replace the corresponding fragments of IgG1 B12, the resulting antibodies being denoted IgG1-CR3022 and IgG1-S309, respectively. It is worth noting that IgG1 is the most abundant IgG subclass in human serum, representing 65% of the total IgG and mediating antibody response against viral pathogens [1]. This molecular engineering process is summarized in Scheme 2. and IgG1-S309 (Fig. 1b), anchored on a gold surface by covalent interaction (Fig. 1c). Details on the construction of the antibodies (Scheme 2) and the sequence of each construct (Figs. S1 and S2) are provided in the Electronic Supporting Information (ESI). A gold nine-layered slab with dimensions 179.5 Â 180.3 Å 2 representing the most stable (1 1 1) facet of the face centered cubic (FCC) gold unit cell [37,38] was built considering 4536 gold atoms per layer and the crystallographic parameter of gold, a = 2.89 Å. The immobilization of the antibodies on gold was performed according to a well-known strategy that considers the activation with carboxylate groups, functionalizing the surface (Scheme 3a) [39,40]. Accordingly, the Lys478 residue of IgG1 B12 which is located at the end of the CH 3 heavy chain H of the F c region, was replaced by the covalent linker displayed in Scheme 3b. Both systems were solvated with 162,597 explicit water molecules and Na + ions were introduced for charge neutralization, the total amount of particles in the simulated systems being 543,221 and 543,223 for IgG1-CR3022 and the IgG1-S309, respectively. Details about the construction of the gold surface, the immobilization of the antibodies and the solvation of the tethered systems are provided in the ESI.

Results and discussion
The constructed models were simulated using AMBER 18 simulation package [41] by applying the force-field parameters described in the ESI. After equilibration and thermalization using the protocols provided in the ESI, 150 ns of NVT conventional MD (cMD) were run for each system. Then, the conformational sampling of the immobilized antibodies was enhanced by performing accelerated MD (aMD) [42], which artificially reduces energy barriers separating different states of the studied system. More specifically, aMD simulations were conducted along 90 ns using an NVT ensemble and starting from 4 snapshots for each studied system (i.e. the 4 initial structures were taken from cMD simulations at 90, 105, 120 and 140 ns). Accordingly, a total of 360 ns of aMD were produced for each tethered engineered antibody. Fig. 2a shows the Root Mean Square Distance (RMSD) along 150 ns of cMD for the two tethered antibodies, which were calcu-  lated with respect to the initial models as constructed using the crystal structures of the CR3022 Fab, the S309 Fab and the IgG1 B12 antibody. The two studied systems were equilibrated after $ 40 ns of production. The averages of RMSD values in the last 110 ns of production take values of 13.5 ± 1.2 Å and 12.8 ± 1.3 Å for IgG1-CR3022 and IgG1-S309, respectively. In both cases, the standard deviation is below 1.5 Å. However, the RMSD values indicate significant distortions with respect to the initial models, which have been attributed to both the fact of combining different crystal structures in the engineered molecules and the restrictions caused by the tethering to the gold surface.
The Root Mean Square Fluctuation (RMSF), which indicates the mobility of the residues during the cMD trajectory, is shown in Fig. 2b. The RMSF fluctuates between 2.45 Å (Cyx984 from the light chain of CR3022 Fab fragment) and 9.95 Å (Asp74 from the heavy chain of CR3022 Fab fragment) for IgG1-CR3022, and between 2.22 Å (Tyr482 from the heavy chain of S309 Fab fragment) and 10.62 Å (Ala923 from the light chain of S309 Fab fragment) for IgG1-S309. A general analysis of the impact of the RMSF in the different regions of the two antibodies is displayed in Fig. 2c, which represents the structure of IgG1-CR3022 and the IgG1-S309 in different colors as function of the RMSF values. As it can be seen, the larger distortions occur in the F c and Fab1 regions for IgG1-CR3022, while they are essentially located in Fab1 for IgG1-S309. The fact that the distortions of the SARS-CoV-2 neutralizing antibodies occur in the Fab1 fragment and such distortions are higher for IgG1-CR3022 than for IgG1-S309 are supported by analyses of the distance between the gold surface and the center of mass of each Fab arm, denoted d Fab1 and d Fab2 in Fig. S3a. Inspection of the temporal evolution of d Fab1 along cMD trajectories (Fig. S3b) indicates that the difference between the average value from the last 110 ns (57 ± 2 and 77 ± 5 Å for IgG1-CR3022 and IgG1-S309, respectively) and the initial value (70 and 72 Å, respectively) is significantly higher for IgG1-CR3022 than for IgG1-S309. Conversely, average d Fab2 values (91 ± 2 and 90 ± 4 Å, respectively) are very similar to the initial values (87 and 92 Å, respectively) for the two antibodies (Fig. S3c).
The effect of the gold surface on the orientation of the engineered antibodies was analyzed considering both cMD and aMD simulations. Firstly, the tilt angle of the engineered antibodies with respect to the constructed surface was examined using the a angle ( Fig. 3a), which was defined by a gold atom of the slab at the top, the a-carbon atom of the residue used to tether the antibody to the surface (i.e. C a in Scheme 3b), and the a-carbon atom of the residue at the top of the F c region (i.e. Pro233 and Pro241 for IgG1-CR3022 and the IgG1-S309, respectively). The evolution of a along the cMD simulation, which is represented in Fig. 3b, shows that stabilization is reached after $ 120 and $ 95 ns for IgG1-CR3022 and the IgG1-S309, respectively. The tilt angle of the initially constructed models (i.e. 72°and 93°for IgG1-CR3022 and the IgG1-S309, respectively) decreases to a = 54°± 1°and 50°± 1°(averaged over the last 30 and 55 ns, respectively).
The potential of mean force (PMF) free energy profiles for the tilt angle calculated from aMD samplings are displayed in Fig. 3cd, which only show, for the sake of clarity, the region with DG 4.2 kcal/mol. As it can be seen, a significant number of local minima are detected within a DG interval of 1.5 kcal/mol for the two engineered antibodies. The global minimum for IgG1-CR3022 appears at a = 54.1°(green diamond in Fig. 3c)), which is very close to the average value obtained by cMD once the trajectory was stabilized (light blue line in Fig. 3c). However, several local minima destabilized by less than 0.1 kcal/mol with respect to the global minimum (i.e. isoenergetic minima) appears at a = 49.5°and 50.7°, suggesting a tilting of around 40°with respect to the ideal orientation (a = 90°). However, due the flatness of the profile in a range of around 10°(i.e. the DG is lower than 1.0 kcal/mol for a values comprised between 48°and 57°), this value can be considered as a rough approximation. On the other hand, the PMF profile obtained for IgG1-S309 (Fig. 3d) shows similar trends, with a global minimum at a = 49.4°(green diamond) and several local minima with DG < 0.1 kcal/mol surrounding it. These results indicate that the tilt angle of anchored IgG1-CR3022 and IgG1-S309 are roughly estimated at around 50°-55°, which are quite similar to that obtained for IgG1 protein at room temperature (a = 66°) [43]. The flexibility and distortion of the Fab1 and Fab2 arms with respect to the Fc region were determined by measuring the b and c angles (Fig. 4a), respectively. These were defined by the C a atom of the residue used to tether the antibody to the surface (Scheme 3b), and the C a atom of the residue on the top of the F c region and the center of mass of the complementary determining regions (CDR), which are the part of the variable chains in IgGs where these molecules bind to their specific antigen, in the corresponding arm. Inspection of the temporal evolution along the cMD simulations (Fig. 4b-c) indicates that the two angles stabilize around 120 and 100 ns for IgG1-CR3022 and IgG1-S309 (green diamonds) deviate from the such average by only 2°for IgG1-CR3022, and by 3°and 1°, respectively, for IgG1-S309, corroborating the cMD results. An important structural feature that deserves discussion is the asymmetry found in the orientation of the two arms, which was defined by the difference between the two angles (D = jbcj). This asymmetry is significantly higher for IgG1-CR3022 than for IgG1-S309. Thus, cMD simulations predicted D = 53°and 8°for IgG1-CR3022 and IgG1-S309, respectively, similar values being also obtained from aMD when the corresponding global minimum are compared (i.e. D = 49°and 10°for IgG1-CR3022 and IgG1-S309, respectively). However, X-ray crystallography [44] and cryo-electron tomography studies in solution [45] reported that the arrangement of Fab arms relative to the Fc is not symmetrical in IgG1 molecules, which was attributed to their flexibility and dynamics. Indeed, literature shows a great variabil- ity in Fab-Fc angles [46]. Detailed analyses of cryo-electron tomography experiments on individual IgG molecules allowed determining the probability distribution of Fab-Fc angles at equilibrium [47]. Results showed that the Fab-Fc angles are approximately uniformly distributed within the two limiting values, 15°and 128°, which reflect the excluded steric constraints imposed by the molecular structure. However, analysis of IgG1 antibodies by individual particle electron tomography suggested that the orientation of the two Fab arms with respect to the Fc tend to be symmetric with D values close to that observed for IgG1-S309 [48].
An important aspect that deserves deeper analysis is the correlation among a, b and c angles, which has been analyzed considering the results from aMD simulations. c. This feature suggests that the movements associated to the flexibility of the Fab arms and the interactions with the gold surface are independent for Fab1 and Fab2. The radius of gyration (R g ) is an indicator of compactness and a measure for flexibility. Analysis of the temporal evolution of the R g for whole proteins indicates that this parameter stabilizes around 20 ns of cMD (Fig. 6a). The values obtained for IgG1-CR3022 and IgG1-S309 by averaging over last 130 ns of cMD trajectories, which are 44.5 ± 0.6 and 48.6 ± 0.6 Å, respectively, reflects that the latter is more similar with the experimental values reported for the IgG1 family (i.e. 49.4 ± 0.4 [49] and 47.9 ± 2.9 Å [50], which were determined by SAXS).
Analyses of the R g values calculated for the Fc, Fab1 and Fab2 fragments separately (Fig. 6b-c) evidence that the difference among the values found for the whole antibodies are mainly due to the Fc, which is more compact for IgG1-CR3022 than for IgG1-S309 (i.e. 26.4 ± 0.6 and 28.1 ± 0.7 Å, respectively, averaged over the last 130 ns). Conversely, the Fab fragments are very similar in all cases, independently of the antibody or the arm orientation (i.e. 24.4-24.5 Å). PMF profiles determined for F c and Fab fragments from aMD simulations, which are displayed in Fig. 6d-f, corroborate the structural features obtained from cMD. Additionally, these results provide some interesting findings that deserve consideration. Firstly, the potential well obtained for the Fc of IgG1-S309 is not only shifted with respect to that of IgG1-CR3022 but also is narrower (Fig. 6d). Similarly, the well found in the PMF of the Fab2 for IgG1-S309 is significantly narrower than those obtained for the Fab1 of the same engineered antibody and for the two Fab arms of IgG1-CR3022. These differences, which have been ascribed to the primary structures of the CR3022 and S309 fragments, explain the higher compactness of IgG1-CR3022. On the other hand, Fig. S4 displays illustrative 2D PMF maps representing the R g of the whole antibody, and of the F c or the Fab fragments against the a, b and c angles. As it can be seen, the R g values are more restricted for IgG1-CR3022 than for IgG1-S309. Furthermore, the restrictions of the former are practically independently of the a, b and c angles.
Other important parameters are those relative to the orientation of the Fab1 and Fab2 branches with respect to the surface, which are described by the u and h dihedral angles. Those dihedrals are defined by a gold atom of the superficial slab, the C a atom of the residue used to tether the antibody to the surface (i.e. the residue replacing Lys478 in IgG1 B12), the C a atom of the residue on the top of the Fc region (i.e. Pro233 and Pro241 for IgG1- CR3022 and the IgG1-S309, respectively), and the center of masses of the CDR of Fab1 (u) or Fab2 (h), as is schematically sketched in Fig. 7a. Analyses of the cMD trajectories indicate that both u and h stabilize after $100 ns (Fig. 7b-c). The average values, which were obtained considering the last 50 ns of cMD runs were u / h = 80°± 6°/ -160°± 11°and 108°± 4°/ -143°± 2°for IgG1-CR3022 and the IgG1-S309, respectively. This result suggests that the relative disposition between the two arms, which has been defined by N = h u, only exhibits a small dependence on the antibody structure (i.e. N = 120°and 109°, respectively).
PMF plots representing the variation of the relative free energy against u or h are depicted in Fig. 7d-g. Although the position of the global minimum found by aMD is consistent with the average values derived from cMD, the deviation between aMD and cMD is higher than those described for a, b and c. However, the regions around the global minimum are very flat and, therefore, such deviations must be considered with caution. In any case, deviations have a relatively small effect on N that increases to $128°and $110°for IgG1-CR3022 and the IgG1-S309, respectively.
Finally, structural changes induced by the anchoring of the two designed antibodies to the gold surface were evaluated by examining the secondary structure of the specific binding regions located at the Fc, Fab1 and Fab2. The secondary structure was determined using the DSSP program [51], which applies a pattern recognition process of hydrogen bond and geometrical features to assign secondary structure to the residues of a protein. The DSSP algorithm was employed considering the snapshots stored during the first 10 ns of the cMD trajectories that, after averaging, allowed to define the initial secondary structure once the conformational tensions associated to the molecular engineering process are elimi-nated. Besides, the converged secondary structure, which includes changes induced by the immobilization on the gold surface, was obtained by averaging the results derived from the application of the DSSP to the snapshots stored during last 75 ns of the runs. Fig. S5 displays the averaged initial and converged secondary structures of the six specific binding regions identified for the Fc fragment of IgG1-CR3022 and IgG1-S309, which have been labeled as Fc_# (where # ranges from 1 to 6). The sequence (primary structure) associated to each of the six Fc_# for each antibody, which is included in Fig. S5, has been labeled using the residue numberings displayed in Figs. S1 and S2. It is worth mentioning that only structural motifs with average populations higher than 20% have been considered as representative and, therefore, are the only included in Fig. S5. As it can be seen, the secondary structure of the six Fc_# domains was dominated by random, bend, a-helix, turn and antiparallel b sheet motifs for the two antibodies. Comparison between the initial and converged secondary structures indicates major conformational rearrangements (i.e. those that imply changes in the populations of the different motifs that add up to more than 30%) for the Fc_2, Fc_3, Fc_5 and Fc_6 domains of IgG1-CR3022 and for the Fc_2, Fc_3 and Fc_6 domains of IgG1-S309. Fig. 8 displays a complete description of the domains that underwent major conformational changes at the Fc fragment for the two designed antibodies. In general, the most important rearrangement corresponds to the organization of the random motifs into well-defined folded secondary structures. An exception to this behavior was the Fc_2 of IgG1-S309, in which the turn and antiparallel b-sheet motifs of the initial structure transformed into random conformations in the converged structure. Another important feature is that, in general, conformational changes are more drastic for IgG1-CR3022 than for IgG1-S309, which is consistent with the RMSD and RMSF tendencies previously discussed (Fig. 2).
Inspection of the averaged initial and converged structures for the six CDRs identified for the Fab1 fragment of each antibody, which have been labeled as Fab1_# (where # ranges from 1 to 6) in Fig. S6, indicates that only the Fab1_2 and Fab1_3 domains of IgG1-CR3022 and the Fab1_5 of IgG1-S309 experience accumulated changes in the populations of the different secondary motifs higher than 30%. However, detailed analysis of the secondary motifs involved in such three domains (Fig. 8) reflects similarities with the specific binding regions of the Fc fragment. Thus, in general the main conformational changes are associated with the conversion of random structures into well-defined regular motifs: antiparallel b-sheet (Fab1_2 and Fab1_3 of IgG1-CR3022) or turn (Fab1_5 of IgG1-S309) structures.
Finally, the secondary structures of the six CDRs of the Fab2 arm, which are labelled as Fab2_# (where # ranges from 1 to 6) in Fig. S7, remained practically unchanged during the MD trajec-tory. Thus, the only domain that experienced conformational changes involving an accumulated population variation higher than 30% was the Fab2_3 of IgG1-CR3022 (Fig. S7). However, detailed analysis of the populations of all identified secondary motifs, which are listed in Fig. 8, indicates that such rearrangements are less important than those observed for Fc and Fab1 domains. More specifically, in the case of Fab2_3, the only relevant conformational change is the transformation of turn and, in less extension, random structures into bend motifs.

Conclusions
In this work we presented classical MD simulations of two antibodies, which were designed to detect SARS-CoV-2, immobilized on a gold (1 1 1) surface through a covalent linker. The IgG1-CR3022 and IgG1-S309 antibodies were engineered by replacing fragments of the IgG1 B12 antibody by other fragments of CR3022 and S309, respectively. Those fragments were found to bind the RBD of SARS-CoV-2, playing a key role in neutralizing the coronavirus. The choice of the gold surface was based on its fre- quent utilization as solid support in immunosensors. Previous works showed that the performance of manufactured antibodybased sensors strongly depends on the influence of the solid support and the orientation and structural stability of the immobilized antibody. Accordingly, in this work cMD and aMD simulations were performed to examine and compare the stability of the engineered antibodies when tethered on a solid gold surface. Although cMD simulations showed that, in average, the tilting was 4°higher for IgG1-S309 than for IgG1-CR3022, the rest of the analyzed structural and geometric parameters underwent higher distortions for the latter than for the former. Consistently, the calculated PMF profiles evidenced, not only that the positions of the minima obtained by aMD were very close to those obtained by averaging cMD values, but also that the potential wells were in general narrower for IgG1-S309 than for IgG1-CR3022. These observations are in agreement with the fact that the RMSD and the RMSF were smaller for IgG1-S309 than for IgG1-CR3022 and with the CDR secondary structure analyses at the beginning and at the end of the simulations. Those analyses showed that changes were more pronounced for the latter than for the former. Overall, these results indicate that between the two engineered antibodies, the IgG1-S309 is the most stable from a structural point of view, preserving the conformation of the neutralizing S309 antibody, and the least affected by the tethering to the gold surface.
The sensitivity and specificity of immunosensors largely depend on the structure of the antibody fragments used to recognize the target pathogen. The structure of such fragments should not be affected by the molecular engineering process used to design the proposed antibody or by the substrate used for the antibody immobilization. Results derived from this study allow to conclude that, although the IgG1-CR3022 and IgG1-S309 antibodies were engineered using identical molecular recognition principles, the latter is more appropriated for the generation of antibody based sensors for COVID-19 diagnosis.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.