Structural dynamics of COVID-19 main protease.

Based on the importance of protease enzymes in functioning some viruses particularly coronaviridae, we have carried out an in silico investigation on the biologically important, yet unmapped phenomenon of activity and internal dynamics of COVID-19 main protease (Mpro) via applying finite-temperature all-atom molecular dynamics simulations. Temperature quench echoes generated by applying two successive cooling signals have therefore been analyzed in terms of the temperature-temperature correlation function of the protease within the harmonic approximation. An exponentially decaying brand of behavior has been found for the calculated echo depth values with increasing time, which has accordingly led to a much small dephasing time of about 150 fs, revealing a significant anharmonicity and therefore an overall structural stiffness for the COVID-19 main protease.

Introduction

Over the last two decades, a large number of coronaviruses have so far been identified in their natural reservoir host (bats) since the outbreak of severe acute respiratory syndrome (SARS) in 2003 [1], [2], [3], [4]. Recently, a new coronavirus identified as COVID-19 virus has also been discovered as the etiological agent responsible for the ongoing 2019–2020 viral pneumonia outbreak beginning in Wuhan, China [5], [6], [7], [8]. Based on the fact that COVID-19 is an emerging situation, very few investigations have so far been carried out on the virus, and no targeted therapeutics or effective treatments are therefore available at present. A possible way to disrupt COVID-19 function, just like some other viruses such as the SARS coronavirus [9], [10], [11], [12], [13], [14], is targeting its main protease (Mpro) as a key enzyme in mediating replication and transcription of the virus [15], [16].

Although the crystal structure of COVID-19 Mpro (PDB ID: 6LU7 [17]) provides a solid basis for investigating the relationship between its structure and function, several questions regarding many dynamic aspects of its enzymatic mechanism such as substrate binding, orientation, catalysis, product release, and the way these processes take place have not yet been answered.

Proteins are soft materials and their biological functions at the most fundamental level are related to physical motions of their constituent atoms [18]. As a result, a detailed investigation on dynamics of COVID-19 Mpro could be necessary for a complete understanding of its function. The present work has therefore been devoted to inquiring into the activity and internal structural dynamics of COVID-19 Mpro within the normal mode picture of proteins. To this end, we have applied all-atom molecular dynamics (MD) simulations, as a widely used tool in studying structure and function of biological entities, and in unraveling the ways they interact with other bio/nano structures [19].

The present research has been organized as follows. Section 2 describes the applied computational setup and details. In Section 3 the obtained results are presented and being discussed. We finally conclude our results and discussion in Section 4.

Materials and method

All-atom MD simulations have been applied to investigate the internal structural dynamics of the COVID-19 Mpro as implemented in the NAMD computer software (version 2.13) [20] along with the VMD program (version 1.9.4) [21] to post-process the outputs. Initial atomic positions of the protease have been taken from the Protein Data Bank (PDB) with entry code 6LU7. The Mpro has been solvated in an orthorhombic water box with volume of 9.2 × 10.7 × 10 nm3 (Fig. 1 ) under periodic boundary conditions with a unit-cell padding of about 2 nm to decouple the weakest periodic interactions possible between the Mpro replicas. The Na and Cl ions have randomly been distributed in order for the whole system to be electrostatically neutral. The July 2018 update of CHARMM36 topology and parameter files [23], [24] have been used as the force fields to calculate the potential energy surface with the standard TIP3P water model [25]. The switching and cutoff distances of about 1.0 and 1.2 nm have also been used for truncating non-bonded van der Waals interactions, respectively. The particle-mesh Ewald (PME) [26] method with grid dimensions of 96 × 108 × 108 has been applied to deal with long-range interactions. The dielectric constant of the system has also been set to . The system has been minimized for 0.4 ns during which the convergences of total energy of the system (E), total energy of the Mpro (E), number of hydrogen bonds, radius of gyration, root mean square deviation (RMSD), and salt bridge distances have accordingly been achieved, as illustrated in Fig. 2 .

Fig. 1

Schematic cartoon diagram of COVID-19 Mpro solvated in an orthorhombic water box–rendered in VMD using Tachyon parallel/multiprocessor ray tracing system [22]. The α-helices, β-strands, and random coils/turns have been shown in pink, yellow, and blue/white, respectively.

Fig. 2

Convergence diagrams of (a) total energy of the system (E), (b) total energy of the Mpro (E), (c) number of hydrogen bonds, (d) radius of gyration, (e) RMSD, and (f) salt bridge distances during minimization–plotted by Gnuplot (version 5.2) [27]. The convergence values are −362797.8 kcal/mol, −3430.353 kcal/mol, 138, 2.228 nm, and 0.8575 Å, respectively, except those of the salt bridges, which are within the range of Å.

Hydrogen bonds are identified within the hydrogen donor-acceptor distance of 3.0 Å with angle cutoff of 20.0∘. The oxygen-nitrogen distance cutoff for salt bridges has also been set to 3.20 Å. The isothermal-isobaric (NPT) simulations have been performed at K and 1.01325 bar using Langevin forces with the Langevin damping constant of about 2.5/ps along with the Nosé-Hoover Langevin piston pressure control. The integration time step is 1.0 fs. The reliability of our setup has also been tested to a great accuracy as described in Appendix A.

Results and discussion

For proteins, internal dynamics including a wide range of timescales from femtosecond bond stretching vibrations to deformations requiring seconds or more, has so far been of particular interest due to the importance of these motions in biological functioning [28], [29], [30]. Among them are picosecond motions, which exhibit a delocalized character that often involve harmonic displacements of a large number of atoms in different parts of the protein [31], [32], [33], [34].

To probe the internal dynamics of the COVID-19 Mpro, we first have equilibrated the system during a NPT simulation at 310 K. We have used degrees of freedom based on the fact that a N-atom protein has internal normal modes after removing the six degrees of freedom associated with overall translation and rotation. The time dependence of the temperature auto-correlation function (C ) is then estimated so that the instantaneous temperature takes the form . From then on, the thermostat is removed and the next MD simulations are carried out in the microcanonical (NVE) ensemble. We indeed approximate the Mpro as a large collection of weakly interacting harmonic oscillators on the picosecond timescale for which only very short trajectories are required. These oscillators (normal modes) are synchronized through two successive quenching processes (cooling signals) as the required perturbations forcing the system to operate in a phase-coherent state in which normal modes oscillate in phase. At the beginning of each quenching process, zero velocities are assigned to all atoms of the Mpro, while their positions are left unchanged. The first velocity reassignment enforces the phase coherence for the oscillators decaying in time with a characteristic dephasing time τ due to frequency dispersion of the normal modes. The second reassignment after the delay time (τ) serves as the probing signal, which through interference with the previously generated coherence of the normal modes, results in a resonance (echo) in the form of a sharp dip in T(t). Such a signal indeed tests the degree of coherence of the Mpro structure via examining both the depth of echo and the instant of time at which it takes place. These two are among the quantitative characteristics of the coherence associated to internal dynamics of COVID-19 Mpro.

Fig. 3 illustrates the normalized temperature auto-correlation function (C ) of COVID-19 Mpro calculated according to the formula

Fig. 3

Normalized temperature-temperature correlation function (red squares) and the exponentially decaying fit function (blue) of the form with a heavy tail, leading to the auto-correlation decay time of about fs obtained via a nonlinear least-square (NLLS) fit using Xmgrace (version 5.1.25) [35].

Fitting an exponentially decaying function of the form to the calculated simulation data via a least-square method results in fs as the auto-correlation decay time of the COVID-19 Mpro.

Fig. 4 shows the two successive quench signals applied to the pre-equilibrated system for ps as well as the related harmonic approximation prediction in a contrasting fashion.

Fig. 4

The first (synchronization) and the second (probing) quench signals applied to the Mpro calculated for ps over the intervals [0, τ] and [τ, 4τ], respectively. The related temperature echo is detected at 2τ in the form of a relatively sharp dip (red circle) in the normalized temperature T(t)/T. Deviation from the harmonic approximation is evident within [τ, 4τ].

The degree of agreement between our echo calculations and the theoretical prediction of the harmonic approximation, namelycould now be seen, where C (t) in (Eq. (1) [36]) is approximated by and because of quenching. Moreover, the average temperature of the system after applying the first and second cooling signals drops to approximately one half (from 0 to τ), and to one third (from τ to 4τ) of the initial equilibrium temperature (T). This is due to the fact that each particle has a potential energy and immediately begins to move again after each quench, converting part of its potential energy back to kinetic energy exactly according to the equipartition theorem. Such a small timescale of order ps is enough for the Mpro to return to equilibrium, as the figure shows. At the echo of the applied probing signal is seen in the form of a relatively sharp dip in the diagram of normalized temperature (T(t)/T). For t > τ, a clear deviation from the equipartition theorem is also observed, indicating that the system is very sensitive to such external perturbations.

To estimate the dephasing time (τ) of the Mpro, we have calculated the temperature echo for 21 values of the delay time ranging from 0.05 to 2.0 ps every 0.05 ps. The echo depth values have accordingly been calculated via subtracting the minimum temperature of the system from the average value after the second quench as shown in Fig. 5 .

Fig. 5

The calculated echo depth values (red triangles) as a function of delay time exhibiting an exponentially decaying behavior with increasing τ. Fitting an exponential of the form [blue] to the simulation data results in τ ≃ 151.638 fs as the dephasing time of the Mpro. The time difference between any two consecutive data points is 0.05 ps.

Fitting to the calculated echo depth values results in a characteristic dephasing time of about τ ≃ 151.638 fs as an inherent property of the Mpro. The smallness of this value could be inferred from the fact that it is indeed around the timescale of hydrogen bond vibrations in proteins, revealing that the COVID-19 Mpro has a high degree of structural stiffness as a physical property.

Deviation from the harmonic approximation [η(τ)] is also evident and is calculated at (the echo point) according toas illustrated in Fig. 6 for different values of τ.

Fig. 6

Deviation of the normalized temperature from those estimated by the harmonic model at the echo point () according to Eq. (2) for different values of the delay time from to 2000 fs. The η exhibits an asymptotic behavior as τ increases. The time difference between any two consecutive data points is 0.05 ps.

It is seen that as τ increases, deviation from theory [Eq. (1)] becomes larger, though according to a limiting behavior, which asymptotically approaches to about 73.

Internal dynamics of the Mpro therefore exhibits a weak harmonicity with a significant deviation from theoretical predictions within the normal-mode description of proteins. The coherence of atomic motions in the COVID-19 Mpro decays on such a small timescale of about τ ≃ 151.638 fs through nonlinear contributions to forces between the constituent atoms. The reason why the calculated echo depth values (Fig. 5), or the behavior of normalized temperature (Fig. 4), are not so exact as those predicted by the harmonic picture reveals that the system has a weak harmonicity, and this is largely due to the existence of anharmonic contributions arising from torsional, electrostatic, and van der Waals interactions that contribute to the potential acting within the COVID-19 Mpro. The temperature of the system at the echo time is weighted the most by those modes that have the most total energy, i.e., those that reach their classical turning points on harmonic potential wells described by the harmonic approximation. As Fig. 5 shows, the echo depth decreases exponentially with τ, indicating that the temperature of the system at the echo point becomes closer to its average value after the second cooling pulse applied. Therefore, one possible explanation for decrease in the echo depth is that the more the value of τ, the less the number of modes that reach their classical turning points at which accordingly reduces the echo depth. On the other hand, that temperature at the echo point does not go completely to zero is due to the small contribution of those modes that have non-zero kinetic energy at . One, however, cannot establish a correlation between decrease in the echo depth on one hand, and increase in the contribution of these modes as τ increases on the other, because at the quench time, their contribution to kinetic energy is exactly zero. Another possible explanation for decrease in the echo depth can then be provided so that there exists a considerable deviation from the harmonic picture and therefore the potential well each particle within the Mpro experiences significantly deviates from being exact parabola as the value of delay time increases.

Conclusions

In summary, we have investigated the activity and internal dynamics of the COVID-19 Mpro within the normal-mode picture of proteins using all-atom MD simulations. We therefore have generated temperature quench echoes by applying two successive quench signals and have then examined them in terms of the temperature auto-correlation function within the harmonic approximation. Analyzing the exponentially decaying trend of the calculated echo depth values with increasing time has led to a much small dephasing time of about 150 fs, exhibiting accordingly a significant anharmonicity and therefore a total structural solidity for the COVID-19 Mpro.

CRediT authorship contribution statement

Ashkan Shekaari: Conceptualization, Data curation, Investigation, Methodology, Validation, Visualization, Writing - review & editing. Mahmoud Jafari: Project administration, Supervision.

Declaration of Competing Interest

We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome. We also confirm that the manuscript has been read and approved by all named authors, and the order of authors listed in the manuscript has been approved by all named authors.