Structural modelling of SARS-CoV-2 alpha variant (B.1.1.7) suggests enhanced furin binding and infectivity.

The B.1.1.7 SARS-CoV-2 strain that has emerged in the UK in early December presents seven mutations and three deletions on S-protein structure that could lead to a more infective strain. The P681H mutation in the "PRRAR" furin cleavage site might affect the binding affinity to furin enzyme and hence its infectivity. Therefore, in this study, various structural bioinformatics approaches were used to model the S-protein structure with the B.1.1.7 variant amino acid substitutions and deletions. In addition to modelling the binding of furin to the cleavage site of the wild-type and the B.1.1.7 variant. Conclusively the B.1.1.7 variant resulted in dynamic stability, conformational changes and variations in binding energies in the S-protein structure, resulting in a more favourable binding of furin enzyme to the SARS-CoV-2 S-protein.

Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is associated with COVID-19 disease that emerged in Wuhan, China, in late 2019 is closely related to the SARS-CoV (2002) and MERS (2012) (Wang et al., 2020a). Despite the proofreading capacity during viral replication, the fast spread of the SARS-CoV-2 virus has led to a high rate of mutations in the viral proteins (Moelling, 2021; Plante et al., 2021). Consequently, resulting in evolved viral variants/strains that are more efficient in host cell entry (Shang et al., 2020) and evasion of the immune system (Wang et al., 2020b). The current SARS-CoV-2 strain (B.1.1.7) that has emerged in the UK (United Kingdom) in early December and extended to other parts of the world has created an alarming situation (Kirby, 2021; O'Toole et al., 2021). The B.1.1.7 strain consists of seven mutations and three residue deletions on the spike protein (S-protein) (Davies et al., 2021; Ostrov, 2021) compared to the initial SARS-CoV-2 strain that was reported in Wuhan, China (Wu et al., 2020), and the D614G strain that originated in Europe (Haddad et al., 2021; Phan, 2020).

The S-protein facilitates the SARS-CoV-2 viral invasion of the host cell by binding to angiotensin-converting enzyme 2 (ACE2), mainly expressed in the lungs, testes, and kidneys (Mohammad et al., 2020b; Wrapp et al., 2020). For the binding process to occur, the S-protein has to be cleaved by proteases. The SARS-CoV-2 S-protein has a unique polybasic “PRRAR” cleavage site, recognized by the furin enzyme (Coutard et al., 2020). Mutations and deletions in the S-protein sequence can alter the structure, affecting its stability and function, further exacerbating SARS-CoV-2 infectivity (Berger and Schaffitzel, 2020; Korber et al., 2020). Increased SARS-CoV-2 infectivity was observed with the S-protein-D614G mutation initially found in the European population (Phan, 2020; Zhang et al., 2020). Patients infected with SARS-CoV-2 D614G mutation manifested higher viral RNA levels and demonstrated higher titers with pseudoviruses in vitro infection assays, indicating that the D614G mutation is more permissive than the wild type (WT) virus first discovered in Wuhan, China(Korber et al., 2020). Furthermore, the structural analysis showed that the amino acid substitution aspartic acid (D) to glycine (G) at position 614 altered the structure of the S-protein, making it easy for the furin to be cleaved (Eaaswarkhanth et al., 2020; Jackson et al., 2021; Mohammad et al., 2020a; Zhang et al., 2020).

The B.1.1.7 SARS-CoV-2 variant strain exhibits missense mutations (N501Y, A570D, P681H, D614G, T716I, S982A, D1118H) and three deletions in residues H69, V70, Y144 (Davies et al., 2021). Notably, the N501Y mutation is located in the receptor-binding domain (RBD), which interacts with ACE2 and the deletions present on the N-terminal domain (NTD) of the S-protein (Khan et al., 2021b). One of the notable mutations in the B.1.1.7 variant is the P681H, which is a substitution of proline (P) to histidine (H) on the S-protein “PRRAR” at furin cleavage site. Since furin proteolysis occurs at a specific multi-basic sequence (R-[X]-(R/K)-R) (Tian et al., 2011), the substitution of P681H on the furin cleavage site may affect the structure of S-protein, hence enhancing the infectivity of the new SARS-CoV-2 variant. Therefore, it would be of interest to ascertain if all the mutations and deletions present in the B.1.1.7 variant affect the binding to furin or just the P681H mutation.

Recently it was shown that the D614G (European origin) substitution influenced the S-protein stability and enhanced binding of furin to S-protein “PRRAR” cleavage site as compared to the WT discovered in Wuhan, China (Eaaswarkhanth et al., 2020; Mohammad et al., 2020). Hence, in this study, we compared the S-protein-B.1.1.7 variant structure in complex with furin with the S-protein P681H-furin complex to determine if the B.1.1.7 variant or the P681H mutation affected furin cleavage. Therefore, we utilized structural bioinformatics tools to decipher the binding differences and structural-dynamic variations of S-protein-P681H and the B.1.1.7 variant. Furthermore, to understand the stability of the S-protein P681H mutation, compared the S-protein B.1.1.7 variant upon interaction with furin, both states were subjected to 100ns molecular dynamic (MD) simulations. In addition, we calculated the binding affinity of furin to S-protein-P681H and S-protein-B.1.1.7 variant using Generalized-Born surface area molecular mechanics (MM/GBSA). We show that the B.1.1.7 variant allows a tighter binding of the S-protein with furin, which justifies the higher infectivity than the WT strain of SARS-CoV-2.

Methods

Structural stability analysis

The structures of spike protein (PDB ID: 6VSB) (Wrapp et al., 2020) and furin (PDB ID: 4Z2A) (Pearce et al., 2019) were used as models for structural analysis in this report. SWISS-Model (Waterhouse et al., 2018) was used to model the missing amino acids that were not visible in the Cryo-EM structure of the S-protein for both WT-P681H and the B.1.1.7 variant. The DynaMut (Rodrigues et al., 2018) webserver was used to introduce mutations on the S-protein structure.

Protein-protein docking

Protein-protein complexes of WT-P681H and B.1.1.7 variant S-protein (PDB ID: 6VSB) with furin (PDB ID: 4Z2A) were modelled by employing a protein-protein docking approach using the HDOCK server (Yan et al., 2020; Yan et al., 2017), which is based on a hybrid algorithm of template-based modelling and ab initio free docking. Restraint docking was performed by defining the critical residues 682, 683, 684, and 685 of the RRAR site as binding site residues of the S-protein structure (PDB ID: 6VSB). Model 1, with the lowest docking energy score and the highest ligand RMSD (Xue et al., 2016).

MD simulations

Solvated complexes of WT-P681H and B.1.1.7 neutralized by counter ions were subjected to all atoms molecular dynamics simulation using force field (ff18) in AMBER20 simulation program (Salomon-Ferrer et al., 2013; Salomon‐Ferrer et al., 2013). Two-steps each 12000 and 6000 conjugate gradient energy minimization cycles were completed to relax the complexes and remove the bad clashes. The heating each complex was heated with default parameters set as 300 K for 200 ps. For density equilibration, using weak restraint for 2ns at constant pressure was executed. GPU accelerated, 100ns MD for each complex using constant pressure. To control the temperature Langevin thermostat with 1 atm pressure and 300K employed (Zwanzig, 1973), while Particle Mesh Ewald (PME) algorithm to evaluate long-range interactions, respectively, keeping the cutoff distances 10 Å. To treat the covalent interactions involving hydrogens SHAKE algorithm was used (Ryckaert et al., 1977).

Structural-dynamics features were estimated to understand the impact of S-protein-P681H and B.1.1.7 variant on the stability, flexibility, compactness of the protein, motion using CPPTRAJ and PTRAJ(Roe and Cheatham III, 2013). Thermodynamic stability of each complex was calculated as the root mean square deviation (RMSD), while the residual flexibility was evaluated as root mean square fluctuation (RMSF).

Binding energy differences estimation

To compute the binding differences elicited on the protein structures of S-protein-P681H and B.1.1.7, structural frames from the conformational dynamics were used to estimate binding free energy. MM/GBSA methods were employed to calculate the contributed Van der Waal, electrostatic and total binding energies of each complex (Hou et al., 2011). These methods are widely used and reported to have a strong correlation with experimental results (Khan et al., 2020; Khan et al., 2021a; Ylilauri and Pentikäinen, 2013). Each energy term such as vdW, electrostatic, GB and SA was calculated as a part of the total binding energy.

For Free Energy calculation, the following equation was used:

Each component of the total free energy was estimated using the following equation:

Gbond, Gele, and GvdW denote bonded, electrostatic, and van der Waals interactions, respectively. G-pol and Gnpol are polar and nonpolar solvated free energies. The Gpol and Gnpol are calculated by the generalized Born (GB) implicit solvent method with the solvent-accessible surface area SASA term.

Results and discussion

S-protein P681H and B.1.1.7 mutations

Genome sequencing of the new B.1.1.7 SARS-CoV-2 strain presented missense mutations (N501Y, A570D, P681H, D614G, T716I, S982A, D1118H) and three deletions in residues H69, V70, Y144. One of the notable mutations in the B.1.1.7 variant is the P681H, which is a substitution of proline (P) to histidine (H) on the S-protein “PRRAR” furin cleavage site. Since the substitution of P681H on the furin cleavage site may affect S-protein structure, enhancing the infectivity of the new SARS-CoV-2 variant. Here, we compared the S-protein-P681H structure with the S-protein-B.1.1.7 variant structure (Fig. 1 ).

Fig. 1

(A) SARS-CoV-2 S-protein trimeric structure PDB:6VSB (Blue). (B) Overlay of the modeled P681H S-protein (orange) and S-protein 6VSB (blue). (C) Overlay of the S-protein 6VSB (Blue) and modeled B.1.1.7 S-protein variant (teal).

The S-protein-P681H and S-protein-B.1.1.7 variant were modeled with the Swiss Model server (Waterhouse et al., 2018) using the 3.5-Å-resolution cryo-EM structure of the SARS-CoV-2 S-protein (PDB ID: 6VSB) (Wrapp et al., 2020) as a scaffold (Fig. 1A). The S-protein-P681H (Fig. 1B) structure showed a root-mean-square deviation (RMSD) of 0.25 Å between modeled and 6VSB structure. Whereby, the modeled B.1.1.7 (Fig. 1C) variant structure presented a 1.093Å RMSD difference with S-protein 6VSB, indicating that the mutations and deletions have induced substantial conformational changes to the S-protein structure, affecting the function of the S-protein.

Furin docking to P681H and B.1.1.7

Furin-S-protein docking was conducted with HDOCK (Yan et al., 2020; Yan et al., 2017) to elucidate how the conformational changes of P681H and B.1.1.7 S-protein structures affect the binding to furin. Model 1 was chosen (Dudenhoeffer et al., 2019) as the best fit for structural analysis. The S-protein cleavage site is positioned in a surface-exposed loop region, making it accessible for furin binding. Structural analysis of the B.1.1.7 loop region demonstrated a conformational change compared to WT S-protein (Fig. 2 A). The introduction of conformational changes near the binding site can influence the furin cleavage process (Yang et al., 2014). Moreover, surface accessible surface area (SASA) analysis showed that the B.1.1.7 (P681H) HRRAR, in the furin cleavage site residues, demonstrated a SASA at the H (163.075), A (176.744), A (216.679), A (100.593), A (84.162), in contrast, the WT PRRAR furin site showed a SASA of P (137.106), A (119.285), A (115.560), A (7.003), A (71.406). Consequently, the higher exposure of the B.1.1.7 S-protein furin cleavage site, making it more accessible to furin cleavage.

Fig. 2

(A) Overlay of the P681H S-protein (orange) and B.1.1.7 variant (teal), residue mutations are depicted in magenta and deleted residues in blue. (B and B) Modeled furin enzyme (pink) bound to B.1.1.7 (teal) and P6816 (orange) S-protein furin cleavage site.

Molecular dynamic simulations of P681H and B.1.1.7 S-protein bound to furin

To characterize the structural stability and dynamic features at the atomic and spatial resolution, S-protein P681H-Furin and B.1.1.7–furin complexes, we ran 100 ns MD simulations using AMBER 20 package. The impact of S-protein-B.1.1.7 and P681H on the interaction with furin was examined by calculating the Cα-atoms root mean square deviation (RMSD) and root-mean-square fluctuation (RMSF) trajectories.

The stability and dynamics of P681H–furin and B.1.1.7–furin S-protein complexes are presented by RMSD values (Fig. 3 A and B). The P681H–furin complex demonstrated a deviation from 0 to 3 Å in the initial 20 ns. The RMSD, then decreased back to 2 Å at 24 ns, deviating to 4 Å during the simulation from 24-40 ns. Subsequently, from 40 to 100 ns, the S-protein-P681H-furin complex displayed high dynamic fluctuations presenting high structural perturbation, with the RMSD decreasing to 2 Å between 40-60 ns and then increased back to 3 Å at the end of the simulation. On the other hand, the S-protein-B.1.1.7-furin complex (Fig. 3B) exhibited a more stable system, where the complex converged to 2 Å at 5 ns, compared to P681H.Subsequentlythe complex stabilised at 2 Å from 5 to 30 ns. After which, the complex converged to 2.5 Å from 30 to 35 ns, where the B.1.1.7-furin complex fluctuated between 2.5 to 3 Å, for the remainder of the simulation (100 ns) indicating a more stable complex than P681H-furin complex.

Fig. 3

Molecular dynamic simulations RMSD plots at 100 ns simulations of S-protein (A) P681H-furin (orange) and (B) S-protein B.1.1.7-furin (teal).

RMSD results for the B.1.1.7–furin complex revealed a stable system that affected the S-protein stability and, therefore, the binding to furin. For instance the global dynamics of the new variants revealed that mutations which increase the stability are strongly corelated with higher infectivity (Starr et al., 2020). This trend was also reported by a recent study on the B.1.617 variant, where the enhanced infectivity was correlated with the increased infection(Khan et al.).The dynamic structure, as a result of the seven mutations and three deletions present on the S-protein structure, has affected the conformational dynamics of the S-protein and consequently conformational optimization space interaction with furin. Consequently, the B.1.1.7-Furin complex remained stable, thus following the global stability trend and enhances the infectivity.

The RMSF values from the 100 ns simulations of S-protein-P681H and B.1.1.7 bound to furin are depicted in Fig. 4 . The P681H and B.1.1.7-furin complexes demonstrated a similar RMSF profile with the P681H-complex demonstrating higher residual fluctuations neat the loop region between residues 800 and 1000 where furin binds to S-protein. The mutations and deletions on the S-protein B.1.1.7 complex have resulted in a conformational change allowing a more stable complex with furin, as shown from the stable internal dynamics of the B.1.1.7-furin complex. In general, conformational changes that result in higher mobility of the structure might result in a more stable complex when interacting with an enzyme (Jernigan et al., 2020; Ruvinsky et al., 2012). The S-protein–furin binding interaction occurs through conformational selection or an induced-fit mechanism (Tsai et al., 2001, Yang et al., 2014), resulting in a more stable structure between B.1.1.7 and furin than the P681H-furin complex.

Fig. 4

RMSF plot of S-protein P681H-furin (orange) and S-protein B.1.1.7-Furin complex (teal).

Binding energies of furin to P681H and B.1.1.7

The total binding free energy of the S-protein-furin complexes was calculated using the MM/GBSA method (Kollman et al., 2000) (Table 1 ). The binding affinity of furin to B.1.1.7 S-protein presented binding free energy of −57.08±0.73 (kcal/mol) compared with −12.4±0.45 kcal/mol for furin bound to P681H S-protein, indicating that the mutations and deletions present in the B.1.1.7 variant S-protein structure resulted in a tighter binding to furin than the S-protein-P681H mutation.

Table 1

The MMGBSA binding free energy variables (Kcal/mol) of S-protein WT, P681H and B.1.1.7 variant in complex with furin.

Complexes	vdW	Electrostatic	SASA	Total Binding Energy (ΔG)
WT-Furin*	105.35±0.15	559.47±1.28	11.97±0.08	-56.78±0.38
G614-Furin*	-103.67±0.11	-456.58±2.28	-11.65±0.07	-61.90±0.32
P681H-Furin	-62.81±1.13	249.30±4.16	-8.92±0.11	-12.39±0.45
B.1.1.7-Furin	-110.29±0.63	-339.57±5.75	-15.56±0.05	-57.08±0.73

Data presented is taken from previously published results (Mohammad et al., 2020a).

In an earlier study, on the WT-S-protein-furin complex demonstrated a binding free energy of -56.78±0.38 (kcal/mol), whereas the S-protein-G614 showed binding energy of -61.90±0.32 (kcal/mol) to furin (Mohammad et al., 2020a) (Table 1). In the current study, the B.1.1.7-furin complex demonstrated binding energy of -57.08±0.78 (kcal/mol), which is tighter than the WT-furin complex by ~1 (Kcal/mol). However, the B.1.17-furin complex is weaker than the G614-furin complex by ~ 3 (Kcal/mol). As such, this would suggest that furin interaction to B.1.1.7 variant may result in more efficient cleavage of the S-protein and subsequent enhanced viral entry, resulting in a more virulent strain than WT but not the G614 strain. Nevertheless, the B.1.1.7 variant possess multiple mutations, especially on the RBD domain that interacts with ACE2. Whereby, a recent structural study involving B.1.17 variant suggested the N501Y variant in the RBD domain enhanced the binding to ACE2, resulting in a more infective virus strain (Khan et al., 2021b, Khan et al., 2021).

Conclusion

Here in this study, we used biomolecular docking and simulation approaches to compare the binding affinities of furin to S-protein P681H and B.1.1.7 variant structures. The B.1.1.7 variant demonstrated seven mutations and three deletions, with mutation P681H on the furin cleavage site. Henceforth, we analyzed the S-protein-B.1.1.7 structure in complex with furin and compared it with the S-protein P681H complex to ascertain if the variant affected the furin binding P681H mutation affected furin cleavage using computational modelling and structural-dynamics approaches. S-protein- B.1.1.7 variant resulted in a conformational change in the furin cleavage site, subsequently affecting the binding affinity of furin to the S-protein. As a result, the S-protein B.1.1.7 variant may enhance viral entry to the host cell causing higher infectivity. A larger-scale investigation will be critical to build a deep understanding of the furin-S-protein interactions and their functional consequences on the spread of the new B.1.1.7 SARS-CoV-2 variant.

CRediT authorship contribution statement

Anwar Mohammad: Conceptualization, Formal analysis, Methodology, Writing – original draft. Jehad Abubaker: Formal analysis, Writing – original draft, Writing – review & editing. Fahd Al-Mulla: Writing – review & editing, Conceptualization.

Declaration of Competing Interest

All authors have no conflict of interest to declare.