Ca2+-dependent mechanism of membrane insertion and destabilization by the SARS-CoV-2 fusion peptide.

Cell penetration after recognition of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) virus by the ACE2 receptor and the fusion of its viral envelope membrane with cellular membranes are the early steps of infectivity. A region of the Spike protein of the virus, identified as the "fusion peptide" (FP), is liberated at its N-terminal site by a specific cleavage occurring in concert with the interaction of the receptor-binding domain of the Spike. Studies have shown that penetration is enhanced by the required binding of Ca2+ ions to the FPs of coronaviruses, but the mechanisms of membrane insertion and destabilization remain unclear. We have predicted the preferred positions of Ca2+ binding to the SARS-CoV-2-FP, the role of Ca2+ ions in mediating peptide-membrane interactions, the preferred mode of insertion of the Ca2+-bound SARS-CoV-2-FP, and consequent effects on the lipid bilayer from extensive atomistic molecular dynamics simulations and trajectory analyses. In a systematic sampling of the interactions of the Ca2+-bound peptide models with lipid membranes, SARS-CoV-2-FP penetrated the bilayer and disrupted its organization only in two modes involving different structural domains. In one, the hydrophobic residues F833/I834 from the middle region of the peptide are inserted. In the other, more prevalent mode, the penetration involves residues L822/F823 from the LLF motif, which is conserved in CoV-2-like viruses, and is achieved by the binding of Ca2+ ions to the D830/D839 and E819/D820 residue pairs. FP penetration is shown to modify the molecular organization in specific areas of the bilayer, and the extent of membrane binding of the SARS-CoV-2 FP is significantly reduced in the absence of Ca2+ ions. These findings provide novel mechanistic insights regarding the role of Ca2+ in mediating SARS-CoV-2 fusion and provide a detailed structural platform to aid the ongoing efforts in rational design of compounds to inhibit SARS-CoV-2 cell entry.

Significance

n class="Species">SARS-CoV-2, the cause of the pan> class="Disease">COVID-19 pandemic, penetrates host cell membranes and uses viral-to-cellular membrane fusion to release its genetic material for replication. Experiments had identified a region termed “fusion peptide” (FP) in the Spike proteins of coronaviruses, as the spearhead in these initial processes, and suggested that Ca2+ is needed to support both functions. Absent structure- and dynamics-based mechanistic information, these FP functions could not be targeted for therapeutic interventions. We describe the development and determination of the missing information from analysis of extensive MD simulation trajectories and propose specific Ca2+-dependent mechanisms of SARS-CoV-2-FP-membrane insertion and destabilization. These results offer a structure-specific platform to aid the ongoing efforts to use this target for the discovery of inhibitors.

Introduction

The n class="Species">coronaviruses (pan> class="Species">CoVs) are envelope viruses that infect human and animal cells via fusion of the viral envelope membrane with cellular membranes (1). Members of the CoV family include the well-known agents of recent pandemics such as severe acute respiratory syndrome CoV (SARS-CoV), the Middle East respiratory syndrome CoV (MERS-CoV) (2), and the SARS-CoV-2 virus responsible for the currently ongoing COVID-19 pandemic (3,4).

The first point of the connection between the n class="Species">SARS-CoV-2 virus anpan>d the pan> class="Species">human cell is at the angiotensin-converting enzyme-2 (ACE2) receptor (1,5, 6, 7), and is now quite well characterized structurally (8, 9, 10, 11, 12, 13, 14, 15, 16, 17). It is mediated by the Spike protein which is a multidomain homotrimer glycoprotein anchored in the viral envelope membrane (5,18). Each monomer of the spike protein consists of two subunits, S1 which includes the receptor-binding domain (RBD) that attaches to the ACE2 receptor (6,7), and S2, which includes a segment termed “fusion peptide” (FP) that is central to the fusion process leading to cell infection by CoVs (5,19). A catalytic cleavage site (the S1/S2 site) is recognized by a furin-like protease that can separate the two component domains of the Spike monomer (20). Host cell recognition of both SARS-CoV and SARS-CoV-2 involves the ACE2 which binds the RBD on S1. The binding of the S1 domain to the host receptor exposes another cleavage site, the S2′ site on the S2 domain, to processing by cell surface enzymes of the host (primarily TMPRSS2) (7). Cleavage of the S2′ site occurs at the N-terminus of the FP segment, and a “tectonic” change in the conformation of the S2 domain repositions the FP regions of each of the three monomers close to one another to form a “spearhead” of the spike, ready to penetrate the cell membrane while still attached to the S1 trimeric bundle (12,21).

FP insertion is thought to result in destabilization of the target cell membrane (22) as an essential step in allowing the entire virus to penetrate into the cell. In the cell, the FP spearheads the fusion of the viral capsule membrane with cell membranes to release the encapsulated genetic material for replication. Thus far, even the rapid pn class="Gene">ace of discovery of structural information about the SARS-CoV-2 (8, 9, 10, 11, 12, 13,23) has not yielded molecular-level data about the vital step of engagement of the viral spike protein with the membrane, much less the structure of its spearhead. Even less is known at the needed molecular level about the mechanism by which the FPs in the spearhead are involved in the fusion of the cell and virus membranes, despite the wealth of information that has recently accumulated from molecular dynamics (MD) simulations on SARS-CoV-2 models (e.g., see Refs (24, 25, 26, 27, 28, 29)).

Here we present results from an atomistic-level MD investigation of the interactions of n class="Species">SARS-CoV-2 FPs with membranpan>es, anpan>d the conpan>sequenpan>ces of these interactionpan>s onpan> membranpan>e properties. Biophysical anpan>d structure/funpan>ctionpan> studies of the pan> class="Species">SARS-CoV-2 FP (30) as well as of the FPs from SARS-CoV (31) and MERS-CoV (32), which are close homologs of SARS-CoV-2, had produced important insights about the functional significance of specific regions of the FP such as the large number of conserved acidic residues in the two consecutive fragments of the FP termed FP1 and FP2 (residues 816–835 and 835–854, respectively, in SARS-CoV-2), and the LLF hydrophobic motif (19,31,33) (L821/L822/F823 in SARS-CoV-2, see Fig. 1). Studies on SARS-CoV-2 FP (30) and SARS-CoV (31) have suggested that together, the F1 and F2 fragments form a bipartite membrane interaction platform, and that acidic residues in both FP1 and FP2 promote membrane binding through their interactions with Ca2+ ions (30,31). These findings about the functional role of Ca2+ are echoed by results obtained for FPs from other related viruses, such as MERS-CoV (32) and Ebola virus (34). Indeed, in the MERS-CoV FP one of these acidic residues, E891 in the N-terminal (FP1) part (corresponding to E819 in SARS-CoV-2 numbering, Fig. 1 A) emerged as crucial for Ca2+ interactions and fusion-related membrane perturbation effects (32). In contrast to the SARS-CoV FP and SARS-CoV-2 FP, however, isothermal titration calorimetry (ITC) experiments with the MERS-CoV FP showed that the latter binds only a single Ca2+ ion (32). Notably, much of the structural and functional inferences obtained from these studies appears to be in excellent agreement with, and recapitulated by, the recently published results obtained from equivalent experiments with the SARS-CoV-2 FP (30).

Figure 1

(A) Sequence alignment of the FP segments from the spike proteins of SARS-CoV-2 (top), SARS-CoV (middle), and MERS-CoV (bottom) CoVs. The sequences are colored according to ClustalX scheme. The conserved acidic residues are shown in purple, and the conserved Cys residues, forming a disulfide bond (C840 and C851 in SARS-CoV-2-FP numbering), are highlight in orange. FP1 and FP2 fragments of the peptide are marked and labeled. (B) Structural model of SARS-CoV-2-FP used in this study. The acidic residues are shown in licorice and labeled. The Cys residues forming the disulfide bond are shown in space fill. To see this figure in color, go online.

(A) Sequence alignment of the FP segments from the n class="Gene">spike proteins of pan> class="Species">SARS-CoV-2 (top), SARS-CoV (middle), and MERS-CoV (bottom) CoVs. The sequences are colored according to ClustalX scheme. The conserved acidic residues are shown in purple, and the conserved Cys residues, forming a disulfide bond (C840 and C851 in SARS-CoV-2-FP numbering), are highlight in orange. FP1 and FP2 fragments of the peptide are marked and labeled. (B) Structural model of SARS-CoV-2-FP used in this study. The acidic residues are shown in licorice and labeled. The Cys residues forming the disulfide bond are shown in space fill. To see this figure in color, go online.

The involvement of the conserved LLF motif in the FP1 fragment of the n class="Species">SARS-CoV-FP as a critical structural elemenpan>t for membranpan>e insertionpan> anpan>d fusogenpan>icity of the peptide was idenpan>tified from mutagenpan>esis studies (19,31,33). Also, perturbationpan> of the pan> class="Chemical">lipid bilayer by Ca2+-dependent FP-membrane interactions was shown with electron spin resonance experiments to affect the organization of the lipid headgroup and backbone atoms of the lipids at the interaction site, but not the central region of the membrane hydrophobic core (30, 31, 32,34). Still, the loci of Ca2+ binding to CoV FPs and their role in any specific (but unknown) modes of FP interaction with the membrane, had remained undetermined. This has hindered the interpretation of any measurable effects of viral interactions with the membrane, and thus any approaches to mitigate infectivity by targeting this essential region.

We addressed these fundamental unknowns about the n class="Species">SARS-CoV-2-FPs anpan>d their structure-based dynpan>amic mechanpan>isms with extenpan>sive atomistic enpan>semble MD simulationpan>s. From anpan>alyses of the various trajectories, we predicted multiple ways in which pan> class="Chemical">Ca2+ ions engage the FP’s acidic residues singly and in pairs. From massive simulations of spontaneous membrane binding of all feasible variants of Ca2+-bound SARS-CoV-2-FP, we found two preferred modes of membrane penetration. In one of these modes, the inserted hydrophobic residues are F833/I834 from the FP2 region of the peptide which, in the full spike context, are closely tethered to the S2 domain. The second more prevalent mode of penetration involves the N-terminal portion of the F1 segment, which is free after the S2′ cleavage. In this mode, the hydrophobic insertion consists of the highly conserved residues L822/F823 (LLF motif) and is achieved by the binding of Ca2+ ions to the D830/D839 pair and the E819/D820 pair of acidic residues. Notably, our extensive control simulations to evaluate the role of Ca2+ binding in the process showed that membrane binding/insertion of the SARS-CoV-2 FP is greatly diminished overall in the absence of Ca2+ ions. Moreover, when bound, the peptide exerted structural changes in the lipid bilayer whose nature and magnitude showed dependence on the local peptide concentration. Together, these findings provide new mechanistic insights and important molecular-level details about the SARS-CoV-2-FP penetration into the membrane and the structural dynamics underlying the key role of Ca2+ binding in this process.

Methods

MD simulations of the SARS-CoV-2-FP in water

All computations were based on the atomistic structure of the FP (residues 816–854) from the full-length n class="Species">SARS-CoV-2 pan> class="Gene">spike protein model described in Ref (27). In this structure, residues Cys840 and Cys851 form a disulfide bond (see Fig. 1). For these studies, the FP segment from the full-length model was isolated and capped with neutral N- and C-termini (ACE and CT3, respectively, in the CHARMM force field nomenclature). Protonation states of all the titratable residues were predicted at pH 7 using Propka 3.1 software (35).

For the atomistic MD simulations of the n class="Species">SARS-CoV-2-FP in n class="Chemical">water, the peptide was embedded in a rectangular solution box and ionized using VMD tools (“Add Solvation Box” and “Add Ions,” respectively) (36). The box of dimensions ∼90 × 80 × 82 Å included an ionic solution mixture of 49 Na+, 51 Cl−, 2 Ca2+ ions, and ∼18,000 water molecules. The total number of atoms in the system was 54,510.

The system was equilibrated with n class="Gene">NAMD versionpan> 2.12 (37) following a multistep protocol during which the backbonpan>e atoms of the pan> class="Species">SARS-CoV-2-FP as well as Ca2+ ions in the solution were first harmonically constrained and subsequently gradually released in four steps (totaling ∼3 ns), changing the restrain force constants kF from 1, to 0.5, to 0.1 kcal/(mol Å2), and 0 kcal/(mol Å2). These simulations implemented all option for rigid bonds, 1 fs (for kF 1, 0.5, and 0.1 kcal/(mol Å2)) or 2 fs (for kF of 0) integration time-step, PME for electrostatics interactions (38), and were carried out in NPT ensemble under isotropic-pressure-coupling conditions, at a temperature of 310 K. The Nose-Hoover Langevin piston algorithm (39) was used to control the target p = 1 atm pressure with the “LangevinPistonPeriod” set to 200 fs and “LangevinPistonDecay” set to 50 fs. The van der Waals interactions were calculated applying a cutoff distance of 12 Å and switching the potential from 10 Å.

After this initial equilibration phase, the velocities of all atoms in the system were reset and ensemble MD runs were initiated with OpenMM version 7.4 (40) during which the system was simulated in 18 independent replicates, each for 640 ns (i.e., cumulative time of ∼11.5 μs). These runs implemented n class="Gene">PME for electrostatic interactionpan>s anpan>d were performed at 310 K temperature unpan>der pan> class="Gene">NVT ensemble. In addition, 4 fs time-step was used, with hydrogen mass repartitioning and with “friction” parameter set to 1.0/picosecond. Additional parameters for these runs included: “EwaldErrorTolerance” 0.0005, “rigidwater” True, and “ConstraintTolerance” 0.000001. The van der Waals interactions were calculated applying a cutoff distance of 12 Å and switching the potential from 10 Å.

Preparation of a pure lipid membrane for atomistic MD simulations

Using the CHARMM-GUI web server (41), we designed a symmetric n class="Chemical">lipid membrane composed of a 3:1:1 mixture of n class="Chemical">POPC (1-palmitoyl-2-oleoyl-glycero-3-phosphocholine), POPG (1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-(1′-rac-glycerol)), and cholesterol. The total number of lipids in the bilayer was 600. The membrane was solvated and ionized (with 150 mM salt concentration), for a total of ∼163,000 atoms.

This system was equilibrated with n class="Gene">NAMD versionpan> 2.12 following the stanpan>dard multistep protocol provided by CHARMM-GUI, anpan>d thenpan> simulated with unpan>biased MD for 25 ns. After this step, three indepenpan>denpan>t replicates were runpan>, each 75 ns lonpan>g. These runpan>s implemenpan>ted all optionpan> for rigid bonpan>ds, 2 fs integrationpan> time-step, pan> class="Gene">PME for electrostatics interactions, and were carried out in NPT ensemble under semi-isotropic-pressure-coupling conditions, at a temperature of 298 K. The Nose-Hoover Langevin piston algorithm was used to control the target p = 1 atm pressure with the “LangevinPistonPeriod” set to 50 fs and “LangevinPistonDecay” set to 25 fs. The van der Waals interactions were calculated with a cutoff distance of 12 Å and switching the potential from 10 Å. In addition, “vdwForceSwitching” was set to yes.

MD simulations of SARS-CoV-2-FP interactions with lipid membranes

Interactions of a single, or multiple, FPs with n class="Chemical">lipid membranpan>es were started by placing 1, 3, 6, 9, or 12 replicas of the pan> class="Chemical">Ca2+-bound SARS-CoV-2-FPs in the proximity of a bilayer composed of 3:1:1 POPC/POPG/Cholesterol that had been pre-equilibrated for 25 ns as described above. Parallel simulations of a single SARS-CoV-2-FP in the absence of Ca2+ ions were carried out with the same pre-equilibrated 3:1:1 POPC/POPG/Cholesterol bilayer. The starting structure for the FP for the latter runs was the FP model from Ref (27) (see Introduction). The SARS-CoV-2-FP and membrane were embedded in a solution box containing 150 mM Na+Cl− salt concentration. The total number of atoms ranged from 167,000 to 217,000, depending on the number of SARS-CoV-2-FPs in the simulation box.

Each system was equilibrated with n class="Gene">NAMD versionpan> 2.12 following the same multistep protocol described above during which the backbonpan>e atoms of the FP as well as the pan> class="Chemical">Ca2+ ions were first harmonically constrained and subsequently gradually released in four steps. After this phase, the velocities of all atoms of the system were reset, and ensemble MD runs were initiated with OpenMM version 7.4. Table 1 lists the number of replicates and simulation times for each construct. These runs implemented PME for electrostatic interactions and were performed at 298 K temperature under NPT ensemble using semi-isotropic-pressure-coupling, with 4 fs time-step, using hydrogen mass repartitioning and with “friction” parameter set to 1.0/picosecond. Additional parameters for these runs included: “EwaldErrorTolerance” 0.0005, “rigidwater” True, and “ConstraintTolerance” 0.000001. The van der Waals interactions were calculated applying a cutoff distance of 12 Å and switching the potential from 10 Å.

Table 1

Listing of SARS-CoV-2-FP/lipid membrane systems studied with all-atom ensemble MD simulations

No. of SARS-CoV-2-FPs	Ca2+ ions	No. of replicates	Simulation time per replicate (μs)	Total simulation time per system (μs)
1	−	42	1	42
1	+	42	1	42
3	+	60	0.7	42
6	+	60	0.7	42
9	+	60	0.7	42
12	+	60	0.7	42

Columns from left to right show the following items: the number of SARS-CoV-2-FPs per system, presence or absence of Ca2+ ions, number of independent replicates for each construct, simulation time per single replicate, and total simulation time per system. The membrane was 3:1:1 mixture of POPC/POPG/cholesterol, with a total of 600 lipids (see Materials and methods).

Listing of n class="Species">SARS-CoV-2-FP/n class="Chemical">lipid membrane systems studied with all-atom ensemble MD simulations

Columns from left to right show the following items: the number of n class="Species">SARS-CoV-2-FPs per system, presenpan>ce or absenpan>ce of pan> class="Chemical">Ca2+ ions, number of independent replicates for each construct, simulation time per single replicate, and total simulation time per system. The membrane was 3:1:1 mixture of POPC/POPG/cholesterol, with a total of 600 lipids (see Materials and methods).

For all simulations we used the latest CHARMM36 force field for proteins and n class="Chemical">lipids (42), as well as the recently revised CHARMM36 force field for ions which includes nonbonded fix (n class="Gene">NBFIX) parameters for Na+ and Ca2+ ions (43).

Calculation of pressure profiles from MD simulations

Lateral pressure profiles were calculated for selected FP-membrane systems (see Results) as well as for the pure n class="Chemical">lipid bilayer system described above. For the latter, the last halves of the trajectories of the three 75-ns-lonpan>g indepenpan>denpan>t replicates were combined into a single trajectory for the pressure anpan>alysis. Similarly, for the FP-membranpan>e complexes three 75-ns-lonpan>g MD trajectory replicas were accumulated in NAMD and the last halves of these trajectories were concatenated and used for analysis. These simulations used the CHARMM36 force field, were initiated from the corresponding last frames of the OpenMM trajectories and were performed with the parameters used in the last step of the multistage equilibration protocol described above.

The pressure profiles were calculated with n class="Gene">NAMD following the protocol described in Ref (44). The conpan>catenpan>ated trajectory for anpan>alysis for each system comprised ∼26,000 frames with output frequenpan>cy of 2 ps. The trajectory was first cenpan>tered onpan> the average z positionpan> of the terminal methyl pan> class="Chemical">carbons of the lipid chains set at (x,y,z) = (0,0,0). The bilayer was then divided into a discrete number of equal-thickness (∼0.8 Å) slabs in the direction normal to the bilayer plane. The kinetic, bonded, and nonbonded contributions to the pressure were obtained separately from the Ewald contributions, and the two were added to produce the lateral components of the pressure tensor in each slab. The normal component was calculated as described earlier (see Supporting material in Ref (44)) with the formula:where L is the length of the simulation box in the z-direction, and represents the tangential component of the pressure tensor. From here, the pressure in each slab was .

Results

The FP of the n class="Species">SARS-CoV-2 n class="Gene">spike protein (SARS-CoV-2-FP) harbors six acidic residues (E819, D820, D830, D839, D843, and D848) that are conserved in many coronaviruses, such as SARS-CoV and MERS-CoV (see Fig. 1). As the binding of Ca2+ ions was found experimentally to modulate the ability of the FPs from these CoVs to insert into lipid membranes, all these residues had to be considered putative Ca2+-binding sites (30, 31, 32,34). To identify the modes of Ca2+ interaction with the SARS-CoV-2-FP and to elucidate the function-related consequences of these interactions, we first used atomistic MD simulations of SARS-CoV-2-FP in aqueous solution to explore the spontaneous binding of Ca2+ and its effect on the conformation of the peptide (see Materials and methods).

Modes of Ca2+ binding to SARS-CoV-2-FP

Analysis of 18 independent 640 ns MD trajectories amounting to >11.5 μs sampling (Materials and methods), captured multiple events of association of n class="Chemical">Ca2+ ionpan>s with the peptide. Various modes of interactionpan>s betweenpan> the peptide anpan>d pan> class="Chemical">Ca2+ ions were observed in these simulations by monitoring: 1) distances between the Ca2+ ions in solution and the side chains of all acidic residues in the peptide; and 2) pairwise distances between the side chains of all acidic residues (stable binding of Ca2+ is expected to bring the side chains of the coordinating residues close to each other).

As shown in Fig. 2 A, we observed simultaneous association of 2 n class="Chemical">Ca2+ ionpan>s with differenpan>t pairs of residues in the FP (see red anpan>d blue rectanpan>gles) in 4 out of 18 trajectories (Replica IDs 6, 9, 14, anpan>d 16). Inpan> 3 out of 18 simulationpan>s (Replica IDs 1, 10, anpan>d 17 in Fig. 2 A), binding of onpan>ly onpan>e pan> class="Chemical">Ca2+ ion per FP was observed, with various pairs of residues involved. In the remaining trajectories we found instances in which a Ca2+ ion was associated with a single acidic residue in the FP (the events marked by gray-striped rectangles in Fig. 2 A). Because, in such cases, the coordination of the bound Ca2+ is not optimal and leads to low affinity/stability, these systems were not considered for further analysis.

Figure 2

(A) Modes of Ca2+ binding to SARS-CoV-2-FP. The table shows acidic residues implicated in Ca2+ binding in 18 independent atomistic MD simulations (each 640 ns in length). In a particular trajectory, residue pairs simultaneously engaging with the bound Ca2+ are denoted by red or blue rectangles, whereas instances of Ca2+ ion associating with a single acidic residue is depicted with gray-stripped rectangle. The trajectories in which simultaneous binding of two Ca2+ ions to different pairs of residues were observed are highlighted in green (Replicates 6, 9, 14, and 16). The simulations in which a single Ca2+ ion was bound to a pair of acidic residues are shown in yellow (Replicates 1, 10, and 17). (B) Structural snapshots of Ca2+-bound SARS-CoV-2-FP from the simulations in (A) highlighted in green and yellow (labeled according to replica ID). In the snapshots, bound Ca2+ ion is depicted as yellow sphere, and the N- and C-terminal parts of the FP (termed FP1 and FP2, respectively) are drawn in tan and cyan colors, correspondingly. The acidic residues engaged with Ca2+ ions are shown in licorice and labeled. To see this figure in color, go online.

(A) Modes of n class="Chemical">Ca2+ binding to pan> class="Species">SARS-CoV-2-FP. The table shows acidic residues implicated in Ca2+ binding in 18 independent atomistic MD simulations (each 640 ns in length). In a particular trajectory, residue pairs simultaneously engaging with the bound Ca2+ are denoted by red or blue rectangles, whereas instances of Ca2+ ion associating with a single acidic residue is depicted with gray-stripped rectangle. The trajectories in which simultaneous binding of two Ca2+ ions to different pairs of residues were observed are highlighted in green (Replicates 6, 9, 14, and 16). The simulations in which a single Ca2+ ion was bound to a pair of acidic residues are shown in yellow (Replicates 1, 10, and 17). (B) Structural snapshots of Ca2+-bound SARS-CoV-2-FP from the simulations in (A) highlighted in green and yellow (labeled according to replica ID). In the snapshots, bound Ca2+ ion is depicted as yellow sphere, and the N- and C-terminal parts of the FP (termed FP1 and FP2, respectively) are drawn in tan and cyan colors, correspondingly. The acidic residues engaged with Ca2+ ions are shown in licorice and labeled. To see this figure in color, go online.

Overall, the binding stoichiometry in the simulations was ∼1.8, in agreement with the experimentally measured value of ∼1.9 for n class="Species">SARS-CoV-2 (30) anpan>d ∼1.7 for pan> class="Species">SARS-CoV (31). Furthermore, studies of SARS-CoV have suggested a model (31) whereby one Ca2+ ion binds to the N-terminal part of the peptide (FP1, residues 801–817, corresponding to 816–835 in SARS-CoV-2), whereas another binds at the C-terminal end (FP2, residues 818–837, corresponding to 836–854 in SARS-CoV-2) (see Fig. 1 A). In this respect, it is interesting to note that among the seven simulations of SARS-CoV-2-FP in which Ca2+ ions engaged with pairs of acidic residues (replicates highlighted in green and yellow in Fig. 2 A), the most common associations by far involved the pairs of residues E819/D820 and D839/D843, each of which occurs three out of seven times. D830 was involved in four out of seven replicas but in different pairs (D830/D839, D820/D830, and D830/D848), and D848 was observed in two out of seven replicas in different pairs (E819/D848 and D830/D848). Structural representations of these different modes of Ca2+ binding are illustrated by the snapshots shown in Fig. 2 B. From the time evolution of the secondary structure along the FP sequence in the concatenated trajectory from the replicas of FP simulations (Fig. S1 A), we see several regions in which the helical structure is sustained for longer times (the N-terminus, and the 835–845 residue stretch) whereas folding/unfolding occurs in others. This is evident in Fig. S1 B showing the fraction of frames in which a particular FP residue is part of a helix. We note that these findings agree well with recently published results from circular dichroism (CD) experiments on SARS-CoV-2-FP showing that the FP is mostly unfolded in solution (30).

Association of the Ca2+-bound SARS-CoV-2-FP with the lipid membrane

We investigated with MD simulations the spontaneous binding to n class="Chemical">lipid bilayers of each of the sevenpan> structural models of the pan> class="Chemical">Ca2+-bound SARS-CoV-2-FP in Fig. 2 B. Each of these structures was placed initially next to a membrane bilayer model composed of 600 lipids in a 3:1:1 mixture of POPC/POPG/Cholesterol that mimics the lipid composition used in recent experimental Ca2+ binding assays in cognate systems (30,34). For each of these seven starting points the MD simulations were run in six independent replicates of 1 μs length, amounting to 42 μs total run time (see Table 1). As described in Materials and methods, a set of separate control MD simulations of the SARS-CoV-2-FP with no Ca2+ ions bound was run with the same lipid membrane under otherwise identical conditions.

Fig. S2A shows the frequency of associations with the n class="Chemical">lipid membranpan>e for each residue of pan> class="Species">SARS-CoV-2-FP in the 42 atomistic MD simulations of the Ca2+-bound peptide. This measure was calculated separately for each group of six MD runs initiated from the same replica (identified by its Replica IDs from Fig. 2 B) and illustrates the preferred modes of peptide-membrane association. A residue was considered to be in contact with the membrane if the z-directional distance between its C atom and the neighboring lipid phosphorus atoms (P-atoms) was <4 Å (see Fig. S2 captions for more details).

The contact frequency results show that the modes of association with the membrane involve different structural segments in the seven replicas that represent the different modes of n class="Chemical">Ca2+ binding shownpan> in Fig. 2 B. Thus, Rep6 anpan>d Rep9 enpan>gage with the bilayer mostly via two segmenpan>ts: the middle regionpan> of the peptide (residues 830–840) anpan>d the C-terminal part (see also Fig. S3). Rep10 and Rep14 show associations through the middle region (residue 825–835) and the N-terminal segment, respectively. Rep16 engages with the membrane through multiple regions (N-terminal, middle, and C-terminal). Lastly, Rep1 and Rep17 exhibited relatively weak associations with the bilayer.

Interestingly, in all cases, the peptide segments interacting with the membrane contained residues that were coordinating n class="Chemical">Ca2+ ionpan>s (Fig. S3) suggesting that the pan> class="Chemical">Ca2+ ions enhance the interactions. In turn, MD simulations of the SARS-CoV-2-FP-membrane system without bound Ca2+ ions registered only sporadic and much weaker binding of SARS-CoV-2-FP to the membrane (Fig. S4).

To obtain a measure of the timescales of the interactions characterized in Fig. S2, we first converted all trajectories into time-series of bound and unbound states (Fig. S5, A–G) and then calculated lifetimes for each bound state. The distribution of the lifetimes (Fig. S5 H) identified relatively long-lived complexes, ∼850 and ∼770 ns, respectively in trajectories from Rep9 and n class="Gene">Rep14 (marked by red and orange stars in Fig. S5). To understand the reason for the observed long lifetimes in these trajectories, we analyzed in more detail their modes of insertion into the n class="Chemical">lipid bilayer.

Two major modes of membrane penetration by Ca2+-bound SARS-CoV-2-FP

The extent of peptide penetration into the n class="Chemical">lipid bilayer was quanpan>tified using two complemenpan>tary approaches. Inpan> onpan>e, we monpan>itored z-coordinate of the C atom of each residue in the pan> class="Species">SARS-CoV-2-FP and considered a residue inserted into the membrane if the z-distance between its C atom and the second carbon atom in the tail of a POPC lipid (atom C22 in CHARMM36 notation) was <5 Å. Using this protocol, we found that the simulations initiated from Rep9, Rep10, and Rep14 resulted in the highest frequency of insertion (Fig. S2 B), driven by either the middle segment of the peptide (as in Rep10 and Rep14) or the N-terminal segment (in Rep9). The simulations initiated from the other replica models, either did not show insertion (Rep1), or indicated only transient membrane penetration by the C-terminal segment of the FP (Rep6, Rep16, Rep17). The latter mode of insertion is not feasible mechanistically in the context of the SARS-CoV-2 spike protein because only the N-terminal end of the FP is freed by the S2′ cleavage while the C-terminal end of the FP remains part of the spike.

Importantly, our analysis revealed that the trajectories in which the membrane-bound state is most persistent (see Fig. S5, runs marked by red and orange stars) also achieve the largest extent of membrane insertion (Fig. S2 C, trajectories 13 and 29), but in two distinct modes of membrane penetration seen in Fig. 3 A and Fig. S6. In one (Mode 1), the peptide inserts the F833 and I834 pair of hydrophobic residues. In the other (Mode 2), the bilayer is penetrated by a different pair of hydrophobic residues, L822/L823 in the n class="Gene">N-terminal segment of the peptide.

Figure 3

(A) Plots of the calculated probabilities for finding the following: a POPC lipid headgroup nitrogen atom (N), a phosphorus atom (P), a lipid tail C22 carbon (C22), any atom of the SARS-CoV-2-FP (FP), and of the F833-I834 pair, inside 2-Å rectangular slabs along the membrane-normal z axis. The data were generated from analysis of one of the MD trajectories in which F833-I834 residue pair inserted into the lipid bilayer. The snapshot shows the peptide-bound bilayer in the orientation that matches the direction of the x axis of the plots in the panel below. The P-atoms are shown as silver spheres, and the rest of the nonhydrogen lipid atoms are shown as lines. (B) Close-up view of the system from (A), illustrating insertion of F833-I834 residue pair (in licorice and labeled). Color-code of the peptide is the same as in Fig. 2B. The acidic residues engaged with Ca2+ ions (yellow spheres) are rendered in red sticks. (C) Plots of the calculated probabilities listed in (A) but for the L822/F823 pair. The data for the plots were generated as in (A). (D) Close-up view of the system from (C), illustrating the insertion of L822/L823 residue pair (in licorice and labeled). Color-code is the same as in (B). See also Fig. S3. To see this figure in color, go online.

(A) Plots of the calculated probabilities for finding the following: a n class="Chemical">POPC pan> class="Chemical">lipid headgroup nitrogen atom (N), a phosphorus atom (P), a lipid tail C22 carbon (C22), any atom of the SARS-CoV-2-FP (FP), and of the F833-I834 pair, inside 2-Å rectangular slabs along the membrane-normal z axis. The data were generated from analysis of one of the MD trajectories in which F833-I834 residue pair inserted into the lipid bilayer. The snapshot shows the peptide-bound bilayer in the orientation that matches the direction of the x axis of the plots in the panel below. The P-atoms are shown as silver spheres, and the rest of the nonhydrogen lipid atoms are shown as lines. (B) Close-up view of the system from (A), illustrating insertion of F833-I834 residue pair (in licorice and labeled). Color-code of the peptide is the same as in Fig. 2B. The acidic residues engaged with Ca2+ ions (yellow spheres) are rendered in red sticks. (C) Plots of the calculated probabilities listed in (A) but for the L822/F823 pair. The data for the plots were generated as in (A). (D) Close-up view of the system from (C), illustrating the insertion of L822/L823 residue pair (in licorice and labeled). Color-code is the same as in (B). See also Fig. S3. To see this figure in color, go online.

To establish the specific depth of these insertions in relation to n class="Chemical">lipid headgroup anpan>d backbonpan>e regionpan>s, we calculated the number denpan>sities of all atoms of the peptide alonpan>g the membranpan>e-normal z axis (i.e., the probability of finding the peptide atoms within rectanpan>gular slabs alonpan>g the z-directionpan>, see descriptionpan> in Fig. S7 legenpan>d). These probability distributionpan>s were thenpan> compared to those for atoms of the pan> class="Chemical">POPC lipids at different positions, e.g., the headgroup nitrogen and phosphorus atoms (N- and P-atoms, respectively), the lipid tail second carbon atoms (C22-atoms) to identify the depth of penetration relative to the membrane thickness. Fig. S7 shows that in many trajectories the density of protein atoms overlaps with that of the lipid headgroup atoms, consistent with our results in Fig. S2 A that the peptide adsorbs on the lipid bilayer. But, in addition, the analysis shows that in the trajectories of replicas exhibiting deep penetrations and the longest lifetimes of bound states, the protein density overlaps with that of the C22-atoms (trajectories 13 and 29, marked by red and yellow stars in Fig. S7). This indicates that the peptide penetrates relatively deeper into the hydrophobic region of the bilayer.

Investigation of the z-directional density profiles of each residue of the peptide confirmed the deep insertion of the same residues identified from the insertion count analysis described above, namely, F833 and I834 in trajectory 13 and L822/F823 in trajectory 29. This is shown by the overlap of densities of the n class="Chemical">F833-I834 pair (in Fig. 3 A for trajectory 13) anpan>d the L822/F823 pair (in Fig. 3 C for trajectory 29) with those of the pan> class="Chemical">lipid N-, P-, and C22-atoms. Thus, we find that the inserted hydrophobic residues partition between the P and C22 atoms of the membrane, approximately equidistantly (i.e., ∼3–4 Å from each of them) along the membrane-normal z axis. That these pairs of SARS-CoV-2-FP residues partition inside the membrane on the level of the C22-atoms (see also the structural snapshots in Fig. 3, B and D), agrees with results from electron spin resonance measurements of membrane interactions of FPs from SARS-CoV (31), MERS-CoV (32), Ebola virus (34), as well as from the recent findings for SARS-CoV-2 (30). The experimental results showed that the lipid bilayer organization is perturbed in the area of the lipid headgroups and backbone, but that the perturbation does not spread to the middle of the hydrophobic core of the membrane (see below for more details). That the major mode of membrane insertion for the SARS-CoV-2-FP involves the pair of residues L822/F823 is of particular importance because these belong to the highly conserved LLF motif that was shown to be critical for interactions with lipid membranes and for infectivity of SARS-CoV-FP (31).

Interestingly, we found that the L822/F823 membrane insertion of the n class="Chemical">Ca2+-bounpan>d FP was accompanpan>ied by a structural chanpan>ge in the pan> class="Gene">N-terminal region which adopts a beta sheet conformation (Fig. S8). This finding agrees with the findings from the recently published results of CD experiments showing a (4–7%) probability of the N-terminus of the SARS-CoV-2-FP folding into an extended conformation in the presence of lipid vesicles (30).

Modes of membrane penetration relate to Ca2+-binding positions in the SARS-CoV-2-FP

The observation that the two main insertion modes were established from trajectories starting with the n class="Species">SARS-CoV-2-FP stabilized by the binding of two pan> class="Chemical">Ca2+ ions—one to the N-terminal residue pair E819/D820, and the other to residue pairs in the C-terminal part D839/D843 (Rep9), or D830/D839 (Rep14)—prompted a closer look at the relationship between the Ca2+-binding poses in the FPs and the modes of membrane insertion. Because the number of insertion events observed was rather small (Fig. S6), we constructed an enhanced set of simulations to probe the interaction of multiple SARS-CoV-2-FPs with the lipid membrane. Simulations of such constructs consider the locally elevated concentration of the FPs when one or more trimeric spike protein associate with the host cell membrane and have the advantage of enhancing the conformational sampling of the system.

The size of the constructs for these simulations was guided by the experimental conditions in studies in which the cognate n class="Species">SARS-CoV-2-FPs (30) anpan>d pan> class="Species">SARS-CoV-FPs interactions with lipid membranes were investigated for a range of peptide concentrations (31). Constructs were therefore built with 3, 6, 9, and 12 SARS-CoV-2-FPs (the highest concentration used in the experiments, a 2% peptide-to-lipid molar ratio (30,31) corresponds to the 12 FPs construct). To build these constructs, replicas of Ca2+-bound SARS-CoV-2-FP structures were selected randomly from the seven models shown in Fig. 2 B and placed on the same lipid membrane that was used for the single SARS-CoV-2-FP simulations at random positions of a rectangular grid. The starting configuration for each peptide concentration consisted of 10 different starting grids (Sim0 to Sim9) populated as detailed in Fig. S9 B for the 12 peptide construct structures illustrated in Fig. S9 A. Each of these starting conditions was simulated in six replicates as detailed in Table 1 and Materials and methods.

The simulations of multiple (3, 6, 9, 12) n class="Species">SARS-CoV-2-FPs interacting with the membranpan>e revealed instanpan>ces of at most four FPs simultanpan>eously inserting into the pan> class="Chemical">lipid bilayer, but most frequently observed were concurrent insertions of one to two FPs (Fig. S10). Importantly, we found that the mode of membrane penetration by SARS-CoV-2-FP via its N-terminal segment (Mode 2 insertion described above) was the most prevalent in these sets of simulations. This is illustrated in Fig. S11, which presents the fraction of trajectory time in which each residue of the peptide is inserted into the membrane in the combined trajectories from the 3, 6, 9, or 12 FP constructs. As shown in Fig. S11, these simulations captured mostly the same insertion modes observed for the single SARS-CoV-2-FP construct (Mode 1, Mode 2, and the C-terminus segment insertion which is possible only in the truncated FP construct). The highest frequency of insertions in all systems was found to involve the L822/F823 pair of residues (i.e., Mode 2).

To examine the relationship between the observed modes of membrane insertion and the initial structural models in the large number of trajectories from the multiple-PF constructs with randomly selected initial structure, we first quantified membrane insertions separately for those peptides that shared the same initial structure (irrespective of their initial pln class="Gene">acemenpan>t onpan> the grid). The histograms of membranpan>e insertionpan>s for the sevenpan> separate sets of data corresponpan>ding to the sevenpan> initial models are shownpan> in Fig. 4 A. The striking observationpan> from this comparisonpan> is that Mode 2 insertionpan> for Rep14 model appears in all four sets of simulations (3, 6, 9, and 12 FPs) and results exclusively from Rep14 and Rep1 models, regardless of peptide concentration (see also snapshots in Fig. 4, B–C). Notably, out of the seven models in Fig. 2, only Rep14 and Rep1 include a Ca2+ ion bound to the D830/D839 pair of residues in the FP2 segment of the peptide pair). Together, these results suggest that the insertion of the N-terminal segment and the conserved hydrophobic motif residues is associated with a specific pattern of Ca2+ ion binding to the SARS-CoV-2 FP. This specific pattern involves the E819/D820 and D830/D839 pairs of residues that facilitate the prevalent Mode 2 of SARS-CoV-2-FP insertion.

Figure 4

(A) The frequency of membrane insertion calculated for each residue of SARS-CoV-2-FP from the MD simulations of 3, 6, 9, and 12 FPs interacting with 3:1:1 POPC/POPG/cholesterol bilayer. The data shown in each row is an average over each of the six starting models (Rep-values) of the peptide, identified as in Fig. 2B. Membrane insertion was defined as described in the caption of Fig. S2B. The conditions that resulted in the insertion of the N-terminal segment of the peptide are highlighted with colored boxes (Rep 1; Rep 14). (B and C) Representative snapshots of FP insertion in the 3 FP and 12 FP systems are shown in (B) and (C), respectively, illustrating the insertion of the N-terminus segment into the membrane. The FP that penetrates the bilayer is in yellow. The Ca2+ ions bound to this FP are shown as golden spheres. The residues coordinating the Ca2+ ion are rendered in licorice. The L822/L823 pair which inserts into the bilayer is rendered in licorice and is labeled. To see this figure in color, go online.

(A) The frequency of membrane insertion calculated for each residue of n class="Species">SARS-CoV-2-FP from the MD simulationpan>s of 3, 6, 9, anpan>d 12 FPs interacting with 3:1:1 pan> class="Chemical">POPC/POPG/cholesterol bilayer. The data shown in each row is an average over each of the six starting models (Rep-values) of the peptide, identified as in Fig. 2B. Membrane insertion was defined as described in the caption of Fig. S2B. The conditions that resulted in the insertion of the N-terminal segment of the peptide are highlighted with colored boxes (Rep 1; Rep 14). (B and C) Representative snapshots of FP insertion in the 3 FP and 12 FP systems are shown in (B) and (C), respectively, illustrating the insertion of the N-terminus segment into the membrane. The FP that penetrates the bilayer is in yellow. The Ca2+ ions bound to this FP are shown as golden spheres. The residues coordinating the Ca2+ ion are rendered in licorice. The L822/L823 pair which inserts into the bilayer is rendered in licorice and is labeled. To see this figure in color, go online.

Response of the membrane environment to the prevalent mode of SARS-CoV-2-FP insertion

To probe the arrangement of the individual membrane constituents around the peptide, we tracked the number of n class="Chemical">POPC, pan> class="Chemical">POPG, and cholesterol lipids in the bilayer near the FP in the single peptide construct trajectories where Mode 2 of insertion was stable (see Fig. 3, C–D). As shown in Fig. 5 A, the environment within 4 Å of the peptide is enriched in POPG and cholesterol, specifically around the N-terminal 820–827 segment, so that the compositions of these two components locally near the peptide increases from the average 0.2 to ∼0.3, whereas the composition of POPC decreases from 0.6 to ∼0.4 (Fig. 5 B). A stabilizing effect of these changes is suggested by the appearance of an electrostatic association of the negatively charged POPG lipid headgroups with the basic side chain of residue K825, which is juxtaposed to the 819/820 Ca2+-binding site (Fig. 5 C).

Figure 5

(A) Time-evolution (y axis) of the number of contacts (i.e., within 4 Å) per SARS-CoV-2-FP residue (x axis) with POPC (left), POPG (middle), and cholesterol (right) in the simulations of the single FP inserted into the membrane via the LLF motif (see also Fig. 3, C and D). (B) Time-evolution of the fraction of POPC (blue), POPG (red), and cholesterol (yellow) molecules in contact with the peptide in the same MD trajectory as in (A). A lipid was considered in contact with the FP if the 3D distance between any lipid atom and any FP atom was <4 Å. Dotted lines indicate the bulk compositions of these lipids in the mixture (3:1:1 POPC/POPG/Cholesterol). (C) A snapshot of the system (as in Fig. 3D) highlighting POPG lipids (in licorice) within 4 Å of the peptide. Residues L822, F823, and K825 are shown in space fill and are labeled. The two bound Ca2+ ions are shown as yellow spheres, and the residues coordinating them are drawn in licorice. To see this figure in color, go online.

(A) Time-evolution (y axis) of the number of contacts (i.e., within 4 Å) per n class="Species">SARS-CoV-2-FP residue (x axis) with pan> class="Chemical">POPC (left), POPG (middle), and cholesterol (right) in the simulations of the single FP inserted into the membrane via the LLF motif (see also Fig. 3, C and D). (B) Time-evolution of the fraction of POPC (blue), POPG (red), and cholesterol (yellow) molecules in contact with the peptide in the same MD trajectory as in (A). A lipid was considered in contact with the FP if the 3D distance between any lipid atom and any FP atom was <4 Å. Dotted lines indicate the bulk compositions of these lipids in the mixture (3:1:1 POPC/POPG/Cholesterol). (C) A snapshot of the system (as in Fig. 3D) highlighting POPG lipids (in licorice) within 4 Å of the peptide. Residues L822, F823, and K825 are shown in space fill and are labeled. The two bound Ca2+ ions are shown as yellow spheres, and the residues coordinating them are drawn in licorice. To see this figure in color, go online.

The effects of these rearrangements in local composition, and the presence of n class="Chemical">Ca2+ ionpan>s being able to interact with the membranpan>e anpan>d the hydrophobic residues at the differenpan>t levels of the proximal bilayer leaflet, are expected to modify membranpan>e properties locally anpan>d be propagated more distanpan>tly. To evaluate the conpan>sequenpan>ces of the observed increase in the local pan> class="Chemical">cholesterol concentration around the peptide we calculated the ordering effect on the lipids in the vicinity of the insertion in the radial zones defined in Fig. 6 A. The calculated deuterium (SCD) order parameters for POPC and POPG lipid tails as well as orientations of PC and PG headgroups in the different zones are shown in Fig. 6, B–E. For the lipids situated in the zones closest to the insertion, we find higher SCD-values, indicating higher order in their hydrocarbon chains that propagates to ∼30 Å from the insertion (Zones 1 and 2). Interestingly, we also find the PG headgroups (but not PC headgroups) in more vertical, “upright” orientation near the peptide (indicated by lower angle values with respect to membrane-normal shown in panels 6F and 6G). This headgroup ordering effect is likely due to the observed electrostatic interactions between the PG headgroups and K825 as it is confined to the first zone closest to the FP.

Figure 6

(A) Definition of radial zones in the membrane around the inserted SARS-CoV-2-FP. A lipid was considered to be in Zone 1, 2, or 3 if its phosphorus atom was within the lateral distance of 15, 15–30, or 30–45 Å, respectively, from the center of mass of the FP. All other lipids constitute Zone 4. (B)–(E) Deuterium order parameters for carbon atoms in the sn1 (B and C) and sn2 (D and E) hydrocarbon chains of POPC (B and D) and POPG (C and E) lipids calculated in the zones defined in (A). The analysis was carried out on the trajectory of 1 FP-membrane system (+1 FP) that was used also for the analysis of the pressure profile (see text, also Fig. 7). For reference, order parameters of the protein-free system are shown in dashed lines. (F–G) Orientation of the POPC (F) and POPG (G) lipid headgroups (P–N, and P–C12 vectors, respectively) with respect to the bilayer normal calculated in the different zones defined in (A). The analysis was carried out on the same trajectory used in (B–E) and for pressure profile calculations Fig. 7). The error bars on the order parameters were calculated from bootstrapping analysis. To see this figure in color, go online.

(A) Definition of radial zones in the membrane around the inserted n class="Species">SARS-CoV-2-FP. A pan> class="Chemical">lipid was considered to be in Zone 1, 2, or 3 if its phosphorus atom was within the lateral distance of 15, 15–30, or 30–45 Å, respectively, from the center of mass of the FP. All other lipids constitute Zone 4. (B)–(E) Deuterium order parameters for carbon atoms in the sn1 (B and C) and sn2 (D and E) hydrocarbon chains of POPC (B and D) and POPG (C and E) lipids calculated in the zones defined in (A). The analysis was carried out on the trajectory of 1 FP-membrane system (+1 FP) that was used also for the analysis of the pressure profile (see text, also Fig. 7). For reference, order parameters of the protein-free system are shown in dashed lines. (F–G) Orientation of the POPC (F) and POPG (G) lipid headgroups (P–N, and P–C12 vectors, respectively) with respect to the bilayer normal calculated in the different zones defined in (A). The analysis was carried out on the same trajectory used in (B–E) and for pressure profile calculations Fig. 7). The error bars on the order parameters were calculated from bootstrapping analysis. To see this figure in color, go online.

Figure 7

(A and B) The pressure profile as a function of the membrane-normal z coordinate in the simulations of a single FP (+1 FP, (A)) and 12 FPs (+12 FP, (B)). For comparison, the pressure distribution in the peptide-free system is also shown (−FP). (C and D) Plots of the probability of finding POPC lipid headgroup nitrogen atom (N), phosphorus atom (P), hydrocarbon chain C22 carbon (C22), and any atom of the SARS-CoV-2-FPs (FP) in 2-Å rectangular slabs along the membrane-normal z axis calculated from the trajectories from simulations of a single FP (1 FP, (C)) and 12 FPs (12 FP, (D)). To see this figure in color, go online.

We next investigated if these FP-induced local structural perturbations in the membrane induce more global effects on bilayer properties, such as on the lateral pressure distribution (45). Pressure profiles were calculated (see Materials and methods) from the single FP trajectories of Mode 2 insertion (Fig. 3, C and D), and from one of the replicates from 12 FPs construct in which two FPs were stably inserted via Mode 2 (Fig. S12). The calculated lateral pressure profiles as a function of the membrane-normal z-coordinate in these two trajectories are shown in panels Fig. 7, n class="Gene">A and B (red lines). For comparisonpan>, the black lines in these panpan>els show the pressure distributionpan> calculated for the peptide-free membranpan>e system. To aid the interpretationpan> of these plots in the structural conpan>text of the bilayer anpan>d the peptide, Fig. 7, C anpan>d D show the positionpan> denpan>sities alonpan>g the membranpan>e-normal for selected pan> class="Disease">lipid atoms (P, N, and C) and FP atoms.

(n class="Gene">A and B) The pressure profile as a function of the membrane-normal z coordinate in the simulations of a single FP (+1 FP, (A)) and 12 FPs (+12 FP, (B)). For comparison, the pressure distribution in the peptide-free system is also shownpan> (−FP). (C and D) Plots of the probability of finding n class="Chemical">POPC lipid headgroup nitrogen atom (N), phosphorus atom (P), hydrocarbon chain C22 carbon (C22), and any atom of the SARS-CoV-2-FPs (FP) in 2-Å rectangular slabs along the membrane-normal z axis calculated from the trajectories from simulations of a single FP (1 FP, (C)) and 12 FPs (12 FP, (D)). To see this figure in color, go online.

The pressure profile in the pure bilayer system follows the expected trends (46): in the bulk solution, where the system is stress-free, the lateral pressure is near zero. Upon approaching the membrane, the lateral pressure deviates from zero and develops both positive and negative peaks in the n class="Chemical">lipid headgroup anpan>d backbonpan>e regionpan>s (10 Å < |z| < 30 Å in Fig. 7). The negative peaks are reflective of attractive forces at the membranpan>e/pan> class="Chemical">water interface that shield the hydrophobic core of the membrane from solvent exposure and thus reduce the bilayer area, whereas the positive peaks are due to repulsive forces that increase the bilayer area. The additional positive peak in the bilayer center (z = 0 in Fig. 7) stems from repulsive forces that arise from entropic losses in dynamics of the lipid tails (46).

The insertion of a single FP peptide has significant effects on the pressure distribution in the leaflet proximal to the peptide, and the profile becomes asymmetric between the two leaflets. Indeed, as shown in Fig. 7 A, both the positive and negative pressure peaks in the FP-inserted leaflet (z > 0, see also Fig. 7 C) are higher than in the pure bilayer, but the pressure is unaltered at the bilayer midplane and in the opposite leaflet. These effects are consistent with the changes in the orientation of the phosphon class="Chemical">lipid headgroups as well as the trenpan>ds in the acyl chain order parameter profiles of pan> class="Chemical">POPC and POPG (Fig. 6) whereby the ordering effect caused by FP is visible across most of the length of the lipid chains except for the last few carbons situated at the bilayer midplane.

The analysis of pressure profiles in the 12 FP system with two FPs simultaneously inserted in the membrane (Fig. 7 B) reveals even more striking differences from the bare membrane system. All peptides, including the two FPs that interact with the membrane, are in the z > 0 region of the simulation box, which makes the pressure in the bulk solution (z > 30 Å in Fig. 7 B) deviate strongly from zero. Importantly, whereas the pressure profile retains similar asymmetry as in the 1 FP system, with the top leaflet exhibiting larger positive and negative peaks than the bottom one, the two inserted peptides affect the pressure distribution in the entire bilayer.

To understand the difference between the membrane perturbations in the single FP compared with the 12 FP construct, we examined the structural properties of the two leaflets separately in these two systems and compared them to those in the peptide-free bilayer. We found that the average order parameter values for n class="Chemical">POPC were similar in the two leaflets of the 12 FP conpan>struct, but lower thanpan> those in the bulk (pure) membranpan>e system (Fig. S13). This chanpan>ge is accompanpan>ied by a slight increase in the average area per pan> class="Chemical">lipid (Fig. S14) and shows an overall deviation from bulk properties, indicative of the presence of stress in the membrane.

As noted earlier and shown in Fig. 6, the structural perturbations induced by a single FP can propagate up to a distance of ∼30 Å. At the highest protein-to-n class="Chemical">lipid ratio used in experimenpan>ts anpan>d simulated in the 12 FP system, the distanpan>ce betweenpan> anpan>y two peptides bounpan>d to the membranpan>e surfpan> class="Gene">ace is less than 60 Å (half of the length in x- and y-directions of the simulated membrane). Thus, in the presence of multiple bound peptides the perturbation ranges overlap so that perturbed bilayer properties cannot relax to their unperturbed values. This results in the accumulation of stress across the entire bilayer. Although the structural changes detected in the 12 FP simulation (Figs. S13 and S14) are small overall, and consistent with a low bilayer tension (between 1.5 and 2 mN/m), they are expected to increase with the number of peptides binding to the membrane. Taken together, these findings show that interactions of FPs from SARS-CoV-2 (or SARS-CoV, which shares the highly conserved FP sequence) induce measurable perturbations of lipid bilayer organization and properties that are dependent on the local peptide concentration and can affect the areas of the lipid headgroups, backbone, and most of the length of the lipid chains.

Discussion

In the exploration of targets for small molecule inhibitors that can prevent n class="Species">coronaviruses like pan> class="Species">SARS-CoV-2 from entering human cells, efforts had been directed mostly toward disrupting the molecular interaction of the RBD of the spike protein with the ACE2 receptor (see Refs (5,47,48) and citations therein). But with the growing experimental evidence for the involvement of the FP region of the spike protein in the interaction with the cell membrane during the entrance of the virus (22), and role of this region in the fusion of the cell and viral membranes (5), we considered the FPs as another possible target for efforts to identify and/or design effective therapeutics against COVID-19. However, unlike the RBD/ACE2 interactions of the SARS- and MERS-CoV, which have been mapped out from structural data at the atomistic level, the structure-function relations underlying the molecular mechanisms of their FP regions in the fusion-ready state of the spike protein, lacks this level of definition and understanding. The structure of the FP in the context of the spike trimer of SARS viruses has been determined to some extent. Indeed, the highest resolution structure of the SARS-CoV-2 spike in this state (23) offers interesting clues for a role of the FP2 region (referred to as the “FP proximal region”) in preparing the recognition of the RBD by the ACE2 receptor. But there is no information about the structure or dynamic properties of the FP in the penetration and “postfusion” state (23) of the spike. The computational analyses presented in this work have addressed this gap in the information needed to target the FP region of the SARS-CoV-2 virus for inhibition of membrane penetration.

Using the described analyses of the trajectories from swarms of MD simulations of the n class="Species">SARS-CoV-2-FP we predicted various feasible pan> class="Chemical">Ca2+-binding positions (“poses”) (Fig. 2). One of the most prevalent poses that emerged from the analysis (see Fig. 2 B) binds two Ca2+ ions, one in the N-terminal region (FP1) at the E819/D820 pair, and the other in the C-terminal region (FP2), at D830/D839. Remarkably, simulations probing the spontaneous binding of the SARS-CoV-2 FP to lipid membranes (started from all the identified Ca2+-binding poses under conditions representing various concentrations of the FP), revealed that binding and penetration are achieved in only two major modes. Moreover, the predominant mode of interaction with the membrane involves the Ca2+-binding pose we had identified to be one of the most likely ones, with the ion bound at E819/D820 and D830/D839. The FP models with this Ca2+-binding pose insert the N-terminal hydrophobic residues L822 and F823—which are part of the conserved LLF motif shown to be of functional importance in SARS-CoV-FP (31)—into the hydrophobic portion of the membrane. This predominant mode is stabilized by several factors including the electrostatic interactions of 1) the Ca2+ bound to E819/D820 with the lipid headgroups and 2) the anionic POPG lipid headgroups with the basic side chain of K825 juxtaposed to the LLF segment. This stabilization suggests an answer to the question raised by the experimental finding that Ca2+ can enhance the interactions between the peptide and the lipid membrane and stabilize the insertion. Although the participation of the charged lipid headgroups in this stabilization is interesting, more computational studies involving other anionic lipids will probe further these mechanistic inferences. Indeed, the agreement with the experimental results in the same membranes with the published data for the FP of SARS-CoV (31) and MERS (32) is very reassuring, and as shown by the recently deposited results in bioRxiv for SARS-CoV-2-FP-membrane interactions (30), the agreement extends to the COVID-19 agent. Notably, experimental studies using POPC/POPG/Cholesterol or POPC/POPS/Cholesterol mixtures have produced similar results regarding the strength of the FP-membrane interactions and the extent of the FP-induced membrane ordering (30,31), suggesting that the charge of the lipid headgroup and not its identity, is the critical parameter here.

Given its prevalence in our simulations and in light of the established role of the LLF motif in n class="Species">SARS-CoV membranpan>e binding anpan>d infectivity (19,31,33), we propose that the membranpan>e-insertionpan> mode involving the L822/F823 pair is likely to be funpan>ctionpan>ally relevanpan>t in creating the fusionpan>-competenpan>t conpan>ditionpan>s. We note that as this manpan>uscript was being finalized, we became aware of a complemenpan>tary work in bioRxiv (49), in which the conpan>formationpan>s of membranpan>e-bounpan>d apo models of pan> class="Species">SARS-CoV-2 (i.e., in the absence of Ca2+ ions) are thoroughly evaluated with MD using a coarse-grained model of the lipid membrane (HMMM). In agreement with our findings, these studies also suggested a membrane-binding mode that involves the N-terminal region of the peptide, but with poses and modes of penetration that relate to the absence of bound Ca2+. Consonant with the experimentally demonstrated Ca2+-binding requirement for efficient insertion, our findings lead to a proposed mechanism for the role of Ca2+ in modulating FP-membrane interactions. This mechanism describes how the preferred mode of Ca2+ binding to the E819/D820 pair supports the preference of the N-terminal mode of insertion, which puts into play the conserved hydrophobic motif. The calculated membrane perturbation by this mode of insertion, discussed further below, supports the role of this mode of insertion in the preparation of the fusion of viral and cell membranes. Our finding that this mode of interaction is stabilized by the binding of the second Ca2+ to the D830/D839 pair, brings to light the interesting observation that the alternative binding of Ca2+ to another set of acidic residues in the C-terminus, such as D839/D843 pair, inhibits the fusion-competent mode of membrane binding. Given the high degree of sequence similarity and conservation among FPs of other coronaviruses such as SARS-CoV and MERS-CoV, we propose that these insights about structure-function relations for SARS-CoV-2, are generalizable to other corona viruses.

Analyses of the effects of the n class="Species">SARS-CoV-2-FP insertionpan> onpan> the dynpan>amic properties of the membranpan>e patch used in the MD simulationpan> trajectories (with a compositionpan> mimicking the experimenpan>tal conpan>ditionpan>s of cognpan>ate measuremenpan>ts (34)) reveals a perturbationpan> pattern that is conpan>sistenpan>t with results from earlier studies onpan> FPs of pan> class="Species">SARS-CoV2, SARS-CoV, MERS-CoV, and Ebola virus (30, 31, 32,34). We find that FP-binding modifies the molecular organization in the area of the lipid headgroups and backbone atoms, as well as in most of the hydrocarbon chain layer; effects on the middle region of the bilayer hydrophobic core are minimal. Interestingly, we observed an enrichment of negatively charged PG lipids and cholesterol molecules near the inserted peptide. As described above, the electrostatic association of PG headgroups with the FP appears to provide additional stability for the LLF insertion mode, whereas segregation of cholesterol has an ordering effect on the lipids around FP that extends to ∼30 Å from the insertion site (Fig. 6). Given that the viral particle will bring multiple FPs to operate in close proximity to each other near the membrane, simultaneous interactions/penetration by multiple FPs will destabilize the bilayer as indicated by our observations for the 12 FP system (Fig. 7). Such bilayer perturbations may be essential for the fusion of the viral and host cell membranes. Future mechanistic studies, involving more complete molecular models of the spike as well as lipid compositions closer to those of physiological cell membranes will be necessary to address these mechanistic hypotheses and to test functional implications of membrane perturbations.

At this stage, however, the results offer a structure-specific platform to aid the ongoing efforts to design inhibitors of virus cell entry at a new and mechanistically important target. Moreover, given the conservation of both sequence and cell-penetration mechanisms of specific n class="Species">CoVs knpan>ownpan> to be pan> class="Species">human pathogens, and the evidence that calcium binding is both essential and shared by FPs of such corona viruses, the known structural similarities and differences in the key regions of the mechanisms we describe should be very useful in seeking specific modes of blocking the penetration and infectivity of this larger set of CoVs.

Author contributions

G.K. and H.W. designed the study. G.K. carried out MD simulations. G.K., A.P. and M.D. performed analyses of the MD trajectories. All the authors interpreted results. G.K. and H.W. wrote the manuscript. All the authors participated in the reviewing/editing of the manuscript and approved the final draft.