Molecular dynamics of Middle East Respiratory Syndrome Coronavirus (MERS CoV) fusion heptad repeat trimers.

Structural studies related to Middle East Respiratory Syndrome Coronavirus (MERS CoV) infection process are so limited. In this study, molecular dynamics (MD) simulations were carried out to unravel changes in the MERS CoV heptad repeat domains (HRs) and factors affecting fusion state HR stability. Results indicated that HR trimer is more rapidly stabilized, having stable system energy and lower root mean square deviations (RMSDs). While trimers were the predominant active form of CoVs HRs, monomers were also discovered in both of viral and cellular membranes. In order to find the differences between S2 monomer and trimer molecular dynamics, S2 monomer was modelled and subjected to MD simulation. In contrast to S2 trimer, S2 monomer was unstable, having high RMSDs with major drifts above 8 Å. Fluctuation of HR residue positions revealed major changes in the C-terminal of HR2 and the linker coil between HR1 and HR2 in both monomer and trimer. Hydrophobic residues at the a and d positions of HR helices stabilize the whole system, with minimal changes in RMSD. The global distance test and contact area difference scores support instability of MERS CoV S2 monomer. Analysis of HR1-HR2 inter-residue contacts and interaction energy revealed three energy scales along HR helices. Two strong interaction energies were identified at the start of the HR2 helix and at the C-terminal of HR2. The identified critical residues by MD simulation and residues at the a and d positions of HR helix were strong stabilizers of HR recognition.

Introduction

In 2012, a new fatal viral disease causing pneumonia and death was identified in Saudi Arabia (Zaki et al., 2012). The newly emerged virus was termed as Middle East Respiratory Syndrome coronavirus (MERS CoV) (de Groot et al., 2013). The infection range comprises the Arabian Peninsula and several countries worldwide (Banik et al., 2015; Choe et al., 2017). The danger of MERS CoV is aggravated by fatal outbreaks documented in South Korea and China (Seong et al., 2016).

Despite several years of MERS CoV circulation, there are still many secrets of virus replication and fusion with host membranes that need more study. The structural approach to revealing changes in virus substructures can be of unique importance in determining viral structural dynamics. However, few molecular dynamics (MD) simulations have been carried out to investigate MERS CoV structural changes and the dynamical aspects of MERS CoV molecular domains (Alfuwaires et al., 2017). The viral membrane fusion protein is a rational target for drug discovery, as inhibition of the viral membrane fusion function can lead to cessation of the replication cycle (Liu et al., 2004; Yao and Compans, 1996; Vanderlinden et al., 2010). This approach proved good efficiency against several viral infections as HIV (Carravilla and Nieva, 2018), SARS CoV (Liu et al., 2009) and respiratory syncytial virus (Mackman et al., 2015).

Viral membrane fusion can be accomplished by fusion of the virus spike with a host cell receptor target (Bosch et al., 2003). In most enveloped viruses, the spike protein is composed of two cleavable protein domains that can be cleaved by proteases. This property was recorded with SARS CoV, MERS CoV and mouse hepatitis virus (MHV) (Xu et al., 2004). However, they show considerable structural differences including the size, composition of fusion proteins and the sites of protein cleavage (Yuan et al., 2017; Wicht et al., 2014). The CoV spike is composed of two proteins, S1 and S2. There are two consecutive events that occur at the start of cell infection. The first step is virus attachment, in which S1 comes into contact with the host receptor. For MERS CoV, dipeptidyl peptidase-4 (DPP4) is the target for binding with host cells (Wang et al., 2013; Kandeel et al., 2014). Soon after attachment, S1 is cleaved by proteolytic enzymes to expose a highly hydrophobic membrane binding domain of S2 (Kirchdoerfer et al., 2016). S2 is the fusion protein that integrates with the host cell membrane; its integration is followed by fusion of the viral and host cell membranes. In MERS CoV and the highly related SARS CoV, S2 is associated with protein fusion process (Liu et al., 2004; Forni et al., 2015). During fusion, major conformational changes occur in S2, forming a six-helical bundle (6HB) of three-stranded coiled coils (Lu et al., 2014). Each S2 subdomain contains two motifs, heptad repeat domain 1 (HR1) and heptad repeat domain 2 (HR2). HR1 forms a homotrimer exposing three hydrophobic pockets on its surface (Xia et al., 2014). S2 HR domains pass through three conformational changes during viral membrane fusion. The first is pre-fusion state, in which both HR1 and HR2 are not bound together. The second is pre-hairpin intermediate state in which 6HB is formed. HR2 packs into the three major hydrophobic grooves of HR1. The last stage is stable hairpin formation, thus bringing the viral and cell membranes into proximity, forming membrane bilayer and starting viral membrane fusion (Gao et al., 2013). When three HR1 motifs align together, the central core is predominantly composed of hydrophobic residues. A HR domain is composed of tandem repeat motifs of seven residues, named from a to g. Of the seven residues, the first (a) and fourth (d) are predominantly hydrophobic or bulky (Gao et al., 2013). This feature is the main forerunner in coiled coil formation and becomes stabilized by the long hydrophobic interface. Previous reports showed that CoV spike is assembled in the form of trimers (Lu et al., 2014). It was reported that there are many unassembled monomers found in the cells as well as on the virion surface (Delmas and Laude, 1990). Trimers are the accepted form of completing the fusion process. The functional and dynamical aspects of discrete spike monomers in virions are still not well understood. In this work, we carried out a comparison of structural dynamics of S2 monomer and trimer from MERS CoV.

Molecular dynamics is a gold standard in the evaluation of protein structural changes and stability (Alfuwaires et al., 2017; Perilla et al., 2015). Quantitative assessment of the changes in protein structure using MD simulation will help in understanding the global and local changes of protein domains or subdomains and support the future design of suitable compounds to modulate protein function. Classical tools such as root mean square deviation (RMSD) and more recent algorithms using global distance test (GDT_TS) and contact area difference (CAD) scores are used to evaluate and compare different structures (Olechnovic et al., 2013). To date, only a few studies have been carried out to investigate the MD of viral membrane fusion in general, and specific studies for MERS CoV are scarce. In this work, we used MD simulation to reveal changes in MERS CoV HR structure during fusion and factors affecting HR stability. MD simulation, energy system stability, RMSD, hydrogen bonding, contact mapping of inter-residue and inter-HR interactions, GDT_TS, and CAD scores were used to evaluate HR stabilization mechanisms. For this purpose, we simulated the MERS CoV S2 protein in the YASARA structure software followed by comprehensive analysis with YASARA built-in analysis macros and webservers for the calculation of global and local changes in distance and contact change measures.

Materials and methods

MD simulation

In order to assess the changes of S2 monomer and trimer structures, two softwares with distinct force fields were used.

MD simulation using YASARA and AMBER force field

Structures of the MERS CoV HRs were retrieved from the Protein Data Bank. Two structures were used in this study, 4MOD and 4NJL. Both structures are similar in sequence and well aligned except for 6 additional residues at the N-terminal region in 4NJL. The software YASARA Structure (version 14.12.2) was used for all MD simulations by opting the use of AMBER14 as a force field. The simulation cell was allowed to include 20 Å surrounding the protein and filled with water at a density of 0.997 g/ml. Initial energy minimization was carried out under relaxed constraints using steepest descent minimization. Simulations were performed in water at constant pressure with temperature at 298 K. In order to mimic physiological conditions, counter ions were added to neutralize the system; Na or Cl was added in replacement of water to give a total NaCl concentration of 0.9%. pH was maintained at 7.4. The simulation was run at constant pressure and temperature (NPT ensemble). All simulation steps were run by a preinstalled macro (md_runfast.mcr) within the YASARA package. Data were collected every 250 ps.

MD simulation using NAMD and CHARMM force field

A molecular dynamics simulation was performed using the CHARMM force field (MacKerell et al., 1998) (version 27) in NAMD (Kalé et al., 1999) with a non-bonded van der Waals cut-off of 12 Å. The monomer and trimer protein were solvated in a cubic TIP3 water box (20 Å water layer). Sixteen Na+ and 12 Cl− (26 Na+ and 14 Cl−) ions were included in the monomer (trimer) case to neutralize the system. Periodic boundary conditions (Jorgensen et al., 1983), a constant temperature of 298 K (controlled by Langevin temperature piston), the NVT canonical ensemble, and the particle-mesh Ewald summation for long range interactions were used. After a steepest-descent energy minimization to remove atomic overlaps, the systems were equilibrated for 0.5 ns, followed by a 50 ns production run with data collection every 2 ps. All simulations were run with SHAKE (Ryckaert et al., 1977) using a 2 fs time step.

Calculation of inter-residue contacts

The contact between HR1 and HR2 residues before and after MD simulation was calculated by YASARA Contact Analyzer. The range of analysed residues included all amino acids of HR2 (L1259-Y1280). During calculation, two sets of results were collected based on the calculated free energy. At first, all contacts were calculated without energy restrictions; then contacts were reanalysed based on a −1.6 kJ/mol (0.38 kcal/mol) contact energy cut-off (Pande and Rokhsar, 1998).

HR1/HR2 inter-residual hydrogen bonds

The changes in H-bonds before and after MD were analysed for HR monomer and trimer by YASARA. The ranges of analysed residues were I997-Q1031 for HR1 (residues in direct contact with HR2 without the linker region) and L1259-Y1280 for HR2.

Calculation of secondary structure content

The secondary structure contents of HR monomer and trimer were analysed before and after MD simulation using the YASARA secondary structure analysis wizard. Comparisons were made based on the percentages of helix, sheet, turn, and coil content.

Global distance test (GDT_TS)

GDT_TS is a common measure of global changes in protein structure. GDT_TS is used to compare the structure similarities between two proteins with identical sequence. In comparison with RMSD, GDT_TS is more accurate in measuring movement of small fragments and changes in flexible termini (Kufareva and Abagyan, 2012). The structures of MERS CoV S2 monomer or trimer were imported to YASARA Structure. The initial structures and those after MD simulation were superimposed. The Critical Assessment of protein Structure Prediction GDT_TS score was calculated over a distance of 1, 2, 4, or 8 Å by the global distance test implemented in YASARA software.

Contact area difference (CAD) score

The CAD score is an important measure for structural changes, providing a measure of change in the contact area between two structures (Grzybkowska et al., 2016; Abagyan and Totrov, 1997). For this analysis, contact MD simulation files were submitted to the CAD score webserver (Olechnovič and Venclovas, 2014). The analysed structures output included all atoms-all atoms (A-A), all atoms-side chains (A-S), and side chains-side chains (S-S). The differences in contacts between two similar proteins can be quantitatively measured and inspected by colour display. The colour coding for superimposed contacts in the structures before and after simulation were red and green colours. Therefore, the changes in contacts between the structures in both S2 monomer and trimer can be visually assessed. Furthermore, local contact area differences can be assessed by evaluation of changes in colour output from CAD server contacts-area plot, where red and blue colour indicates lower or higher contact area differences, respectively.

Results and discussion

Bioinformatics and computational tools are widely used for understanding the functional and structural aspects of microbial proteins (Kandeel and Altaher, 2017; Kandeel et al., 2009; Alnazawi et al., 2017). MD simulation is a widely used technique for understanding structural protein changes in response to different effectors (Alfuwaires et al., 2017; Moore et al., 1998; Shen et al., 2003; Kandeel and Kitade, 2018). In this study, MD simulation was run in a system comprising monomer or trimer of MERS CoV S2 HR. The stability of each system was evaluated by changes in RMSD as well as changes in the system energy. In order to get maximal precision, the MD simulation results were compared from two different software programs by implementing two different force fields, AMBER14 and CHARMM. All MD simulations showed rapid energy stabilization for both HR monomer and trimer. Fig. 1 shows the changes in RMSD for each structure in relation to time in ps. HR trimer showed rapid stabilization at less than 5 ns, having constant low fluctuations in RMSD and remaining around 3 Å over the entire recorded simulation. In contrast, S2 monomer from two structures was less stable, showing high fluctuations in RMSD with major drifts at 25–30 ns (Fig. 1A). Despite the lower RMSD observed for the monomer in 3MOD structure, it shows high fluctuations in RMSD. This indicates that monomer of S2 bears high flexibility and instability, while trimer constitutes the more or less rigid state of S2. This agrees with the prediction models and resolved structures indicating that S2 of SARS CoV (Deng et al., 2006; Bernini et al., 2004) and of MERS CoV could arrange into trimers (Gao et al., 2013). Additionally, the results from NAMD CHARMM run (Fig. 1B) was highly comparable with YASARA AMBER14, indicating conserved features of trimer stability and monomer dynamic nature. Fig. 1C shows the energy during MD simulation and indicates the stability of the trimer at lower energy level.

Fig. 1

Time dependence of RMSD for MERS CoV HR monomer and trimer. Simulation was run for 50 ns. The trace was based on RMSD of α-carbon atom in PDB structures of S2 monomer or trimer from two different structures using AMBER14 (A) and CHARMM (B) force fields (B). The energy during MD simulation is represented in (C).

The changes in RMSD for every amino acid in MERS CoV HR were estimated for trimer (Fig. 2 A) and monomer (Fig. 2B). The crystal structures of monomer (4MOD and 4NJL) showed more or less similar profiles, albeit with some differences in RMSD (Fig. 2B). In S2 monomer, there was more generalized change in RMSD with clear differences at a) the N and C-termini of the HR complex, b) in the middle of the HR1 helix, and c) in the linker between HR1 and HR2. In contrast, the trimeric structures showed different profiles, with major changes in the linker and C-terminal regions and little or no change in other HR regions (Fig. 2A). In addition, most residues in trimer showed low RMSDs of around 1 Å, with a maximum value at 6.2 Å. A large increase in RMSD values was observed for residues in the range from GLY1250 to ASN1256 (RMSD 3–5 Å). S2 monomer showed more dynamic changes, with a peak RMSD exceeding 10 Å and generalized changes of 2–4 Å along the HR residues. Alignments of pre- and post-MD simulation structures for both monomer and trimer are represented in Fig. 2C. The alignment reveals greater stability for trimer (represented by one chain in the lower panel), compared with more dynamic changes in monomer (upper panel), especially in the middle of HR1, the linker region, and at the protein termini.

Fig. 2

RMSD changes in all residues of MERS CoV S2 trimer (A) or monomer (B) after using YASARA software and AMBER14 force field. The data for each monomer is provided sequentially, each monomer starts at residue Leu996. Alignments of pre- and post-MD simulation structures (C) from HR monomer (upper panel) or trimer (lower panel). The pre-MD structure is provided in blue, while the post-MD structure is provided in brown. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The obtained results from NAMD software and CHARMM force field (Fig. 3 ) were almost similar that estimated by YASARA software. This confirms the finding that residues in monomer are highly mobile either within the linker region or within the backbone of HR1 and HR2. The higher RMSD scale (x-axis) in monomer implies generally higher changes in residues in comparison with trimer.

Fig. 3

RMSD changes in all residues of MERS CoV monomer (A) or trimer (B) after using NAMD software and CHARMM force field. Alignments of pre- and post-MD simulation structures are shown above each RMSD/residue plot. The pre-MD structure is provided in light blue, while the post-MD structure is provided in yellow. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

RMSD is a common measure of global changes in protein structure. However, several concerns and uncertainties in using RMSD have been previously raised (Kufareva and Abagyan, 2012). Of special interest are major dynamical changes at the termini of HR domains with large RMSDs, which might result in misestimation of dynamical changes across the whole system. For more accurate consideration of global changes and accurate inclusion of flexible or terminal highly mobile loops, analysis was also performed using GDT_TS. In agreement with RMSD, GDT_TS revealed the stability of MERS CoV S2 trimer (Table 1 ). The percentages of superimposable residues within 1, 2, 4, or 8 Å in trimer were 2- to 3-fold higher than in monomer. This reflects the more dynamic nature of S2 monomer during MD simulation. The GDT_TS scores for monomer and trimer were 40.6 and 74.2, respectively. Therefore, the greater global changes in S2 monomer are decomposed by trimerization.

Table 1

Global distance test (GDT_TS) of MERS CoV S2 monomer or trimer at different cut-off values.

Structure	Percentage of matched atoms that can be superimposed
MERS CoV S2 monomer
Cut-off (Å)
1	13.2
2	19.7
4	27.7
8	41.6
MERS CoV S2 trimer
Cut-off (Å)
1	40.1
2	68
4	90.8
8	98.1

In addition to global changes in HR, specific residue changes were also investigated. The helical component of HR is composed of several repeats of seven residues. The position of residues in these repeats can be termed a, b, c, d, e, f and g. Of special interest are the residues at positions a and d; which are located almost in the center of the HR, are predominantly bulky and hydrophobic, and share in establishing the hydrophobic core of HR. Position a is represented by residues F1012, F1019, V1026, and L1033, while position d comprises residues F1001, M1008, T1015, V1022, and L1036. MD simulation revealed that residues at a and d positions are most stable, having the least changes in RMSD in comparison with the initial structure. The average RMSD after MD simulation for residues at a and d positions (RMSDad) was found to be smaller than the average of all residues. For HR1, the average RMSDad was 4.37 Å and 0.93 Å for HR monomer and trimer, respectively. Similarly, the RMSDad for HR2 was 3.9 Å for monomer and 1.01 Å for trimer. These values are much lower than the general RMSD averages of 4.98 Å for monomer and 1.49 Å for trimer (Table 2 ).

Table 2

RMSD values (Å) for residues at a and d positions of MERS CoV HR.

HR1			HR2
Residue no.	monomer	trimer	Residue no.	monomer	trimer
PHE 1001	2.44	0.92	LEU 1259	5.33	1.33
MET 1008	3.06	1.15	LEU 1262	5.24	0.94
PHE 1012	3.45	0.63	MET 1266	4.01	1.01
THR 1015	4.29	0.75	LEU 1269	3.78	1.01
PHE 1019	4.92	0.84	VAL 1273	2.49	0.92
VAL 1022	5.58	0.96	LEU 1276	2.52	0.90
VAL 1026	5.31	0.98
LEU 1033	5.19	1.14
LEU 1036	5.16	1.06
average	4.37	0.93	average	1.02	3.90

To determine the key factors governing the stabilization of viral HR, the inter-HR1-HR2 contacts were analysed. Residue-residue contacts were also analysed for their energy contributions to HR stabilization. During residue-residue contact calculations, the contact interaction could be significant if the interaction energy was below −1.26 kJ/mol. For the identification of key residues in contacts between HR1 and HR2, the contact value and number of residues were calculated. In all of the analysed data, there was no positive or repulsive energy. After MD simulation, the total number of contacts was increased for the trimer and to a lesser extent in monomer (Table 3 ). Analysis of every HR residue-residue contact revealed three different levels of interaction energy: a) high interaction energy above 10 kJ/mol, b) medium interaction energy of 4–10 kJ/mol, and c) low interaction energy of 1–4 kJ/mol. The high interaction residue contacts occurred at two positions: first, just in proximity to the N-terminal of HR2, between K1021and Q1023 of HR1 and D1261 and L1262 of HR2; second, distal to the C-terminal of HR2, including the interactions between Q994, K1000, D1282, and E1285 (Fig. 4 A). Parallel to the high interaction residues, several lines of medium interaction energy residues were observed (Fig. 4B). These medium interaction residues were distributed at almost regular intervals starting at the end of the linker between HR1 and HR2 (residues E1039 and L1252) and at residues T1257, L1259, L1269, and D1282. Weak interaction energy contacts fill the gaps between the previously described high and medium interaction contacts. This described profile applies to both monomer and trimer. However, in trimer there was an additional high-energy interaction at the start of the linker region. Therefore, it is suggested to consider the interaction energy of residues during the design of new antiviral membrane fusion agents based on short peptides.

Table 3

The secondary structure content, total surface area and residues contacts of MERS CoV S2 monomer and trimer before or after MD simulation.

	Monomer		Trimer
	Before MD simulation	After MD simulation	Before MD simulation	After MD simulation
Helix	71.3	73.6	67.4	72.1
Sheet	0	0	0	0
Turn	3.2	6.2	5.2	2.1
coil	25.6	20.2	27.4	24.8
Total surface area	9564.72	7910.32	17022.02	17243.12
Residue total contacts	481	482	1741	1761

Fig. 4

Residue contacts between MERS CoV HR1 and HR2 with high (A) or medium (B) interaction energies. The interaction energy was calculated by YASARA software.

CAD scores were used to assess the changes in structures after MD simulation, compared to the initial conformation. Similarities and differences in contact areas were plotted on a colour scale of blue, white, and red corresponding to the range from agreement to difference between structures. The superimposed contact map revealed more changes in contacts for monomer of MERS CoV S2 (Fig. 5 ).

Fig. 5

Superimposed contact maps of MERS CoV S2 contacts. Red colour indicates contacts in the initial structure, green colour indicates contacts in the final structure, and yellow colour indicates common contacts in both the initial and final structures. Higher degree of green dots in monomer in comparison with trimer indicates higher changes in monomer after MD simulation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

In monomer, the residues with major changes in contact were ASP1053, ASP1059, GLU1062, SER1064, ARG1067, and GLY1068. In trimer, major contact changes were observed in ARG1067, GLY1068, I1070, and ASN1111. Analysis of A-A contacts revealed a wider area of contact changes in monomer from GLY1045 to LEU1085, while in trimer more restricted distances were seen from GLN1063 to PHE1073. The CAD score was higher for trimer than monomer (Table 4 ). This agrees with the lower RMSD for trimer and indicates low perturbations in trimer and high perturbations in monomer. The residues showing red spots on CAD superimposition plots also had the highest RMSDs in both monomer and trimer. This indicates the feasibility of using RMSD for evaluation of structural changes. Fig. 6 shows local contact area differences plotted on a colour scale for both MERS CoV S2 monomer and trimer before and after MD. Monomer showed more dispersed red spots, indicating larger changes in contact areas. While A-A analysis shows small areas of contact changes, A-S and S-S determinations show larger contact area changes.

Table 4

Contact area difference score (CAD score) of MERS CoV S2 monomer or trimer.

Structure	CAD score
MERS CoV S2 monomer
A-A	0.69
A-S	0.53
S-S	0.24
MERS CoV S2 trimer
A-A	0.78
A-S	0.66
S-S	0.59

Fig. 6

Colour-coded profiles for contact area changes. The colour scale ranges from 0 (blue) to 1 (red). Blue colour indicates lower contact differences. Red colour indicates higher degree of contact area differences. Monomer showed higher differences indicating their variable and dynamic structure. In contrast, trimer was more or less stable by showing lower differences. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The secondary structure content of HR is shown in Table 3. Helices and coils are the major constituents of HR. After MD simulation, the helix% was increased in both monomer and trimer on the expense of coils. Despite the differences recorded between monomer and trimer during MD simulation, little or no significant change was observed in their secondary structure contents. This suggests that the components of HR retain their full helical or secondary structures even before trimerization.

Conclusion

During viral membrane fusion with the cell membrane, the virus spike S2 protein arranges in a coiled coil with its HR2 domain packed into a deep groove on HR1. By MD simulation, we show that monomer is more dynamic and its residues have more positional fluctuation than in trimer. Furthermore, HR2 recognition by HR1 occurs through three levels of energetic interaction, with high, medium, and low energies distributed in parallel patterns along the HR. The hydrophobic residues at the a and d positions of HR helices have the smallest RMSDs. GDT_TS and CAD scores coincide well with RMSD data, supporting the finding that monomer is unstable and undergoes large fluctuations. Based on these results, the design of peptide analogues could consider the energetic and dynamic aspects of HR1 and HR2 interactions. Since discrete or unassembled monomers are found in the cell and in virions, the noticed flexibility and high dynamics of spike monomers might modulate the virus infection process. Additionally, the stable less dynamic trimer might be required in stabilizing the viral-cell membrane hairpin formation in preparation for fusion of virus and cells.

Conflict of interest

The authors declare no conflict of interest.