Evolution and Pathogenicity of SARS-CoVs: A Microcanonical Analysis of Receptor-Binding Motifs

The recent emergence and rapid dissemination of coronaviruses [1, 2], particularly SARS-CoV-1 and SARS-CoV-2, have profoundly impacted global health, emphasizing the urgent need for a detailed understanding of their molecular mechanisms, especially those involved in viral entry and infection [3, 4, 5]. Both SARS-CoV-1, responsible for the Severe Acute Respiratory Syndrome (SARS) outbreak in 2002-2003, and SARS-CoV-2, which triggered the COVID-19 pandemic starting in late 2019, belong to the Coronaviridae family. These viruses share a large, enveloped, positive-sense single-stranded RNA genome, and their origins are closely linked to zoonotic spillovers, likely from bats to intermediate hosts before reaching humans [6]. Despite their similarities, the evolutionary paths of these viruses have led to significant differences in their biological behavior and impact on human health, making the study of their Spike (S) protein critical for understanding their pathogenicity and transmissibility [7].

The Spike protein is a trimeric transmembrane glycoprotein that mediates the attachment of the virus to the host cell surface and the subsequent fusion of the viral and cellular membranes, which is critical for viral entry. The receptor-binding domain (RBD) of the Spike protein is the region that directly interacts with the angiotensin-converting enzyme 2 (ACE2) receptor on host cells, determining the virus’s ability to infect [8, 9]. Within the RBD lies receptor-binding motif (RBM), a key subdomain responsible for the direct interaction with ACE2. This motif is highly conserved in coronaviruses but exhibits subtle variations between different strains and variants, influencing the binding affinity and, consequently, the virus’s transmissibility and pathogenicity [10].

SARS-CoV-1 and SARS-CoV-2 share a common ancestry but differ significantly in the structural features of their Spike proteins, particularly in the RBM. SARS-CoV-1’s RBM binds to ACE2 with high affinity, which correlates with its ability to cause severe respiratory illness [11, 12]. In contrast, SARS-CoV-2’s RBM, although it shares structural similarities with SARS-CoV-1, has evolved to enhance binding affinity to ACE2, contributing to its higher transmissibility [3, 7, 13, 14]. Notably, early variants of concern of SARS-CoV-2, such as the Beta (P.1) and Gamma (B.1.351) variants, exhibit mutations within the RBM that alter key residues involved in ACE2 binding. These mutations, including the N501Y and E484K, not only enhance the binding affinity to ACE2 [15, 16] but also confer the ability to evade neutralizing antibodies generated by natural infection or vaccination [17, 18].

Moreover, the impact of these mutations extends beyond immune evasion [19, 20]. They also influence the Spike protein’s biophysical properties, such as charge distribution, solubility, and structural stability [21]. For instance, the E484K mutation, found in both the Beta ($\beta$) and Gamma ($\gamma$) variants, introduces a positive charge in the RBM, which can disrupt local electrostatic interactions and modify the overall stability of the protein [22, 23]. Similarly, mutations like N501Y enhance the Spike protein’s structural stability, which may contribute to more efficient viral entry [21]. Therefore, the structural and functional differences in the RBM of those SARS-CoVs underscore the complex relationship between viral evolution, protein dynamics, and pathogenicity.

From a physical perspective the folding behavior of RBMs, crucial for their thermostability and Spike interaction potential, can be rigorously analyzed using multicanonical simulations [24, 25, 26, 27], which enable precise calculations of the density of states (DoS). This approach, rooted in D.H.E. Gross’s theoretical framework [28], facilitates microcanonical analyzes that uncovers critical phenomena such as phase transitions [29], often missed in simulations using the canonical ensemble [30]. RBMs, as protein motifs, qualify as “small systems” under this scenario, where the small number of particles such as the limited amino acid residues in RBMs allows surface effects and fluctuations to significantly impact their behavior [31]. Small systems, such as atomic nucleai [32, 33], nanoscale materials [34, 35, 36], and proteins [31, 37, 38] as the Spike RBM motifs, exhibit unique thermostatistical properties, including energy-dependent equilibrium states in the microcanonical ensemble that lack direct counterparts in the canonical ensemble due to the nonconcavity of the entropy function [30, 31, 34]. This nonconcavity can lead to phenomena like negative heat capacity and metastable or unstable states, which canonical methods often overlook [28]. Thus, microcanonical analysis is crucial for accurately capturing these states and understanding the detailed thermostatistics and functional behavior of RBMs.

In addition to folding dynamics, the solubility of proteins in an aqueous environment plays a crucial role in their biological function [39]. The solubility of these peptides can be described using the Poisson-Boltzmann equation (PBE) [40], which accounts for the electrostatic interactions between the protein and the surrounding solvent. By solving the PBE [41, 42] for each RBM variant, we can gain insights into how changes in the amino acid sequence influence these protein’s domains electrostatic potential, and consequently, its solubility [43]. These mutations are also expected to influence the propensity of Spike to form aggregates or interact with other molecules [44, 45], with potential consequences for viral infectivity and immune recognition [10, 46].

The combination of multicanonical simulations [27, 47], microcanonical analysis [28, 31], and electrostatic modeling [40, 48] provides a powerful methodology for understanding the molecular mechanisms that underlie the differences between SARS-CoV-1 and various SARS-CoV-2 variants. By delving into the detailed thermostatistics and solubility properties of the RBM, we aim to uncover the fundamental principles governing the folding, stability, and function of these critical viral motifs. This knowledge is essential not only for understanding the molecular basis of viral evolution and pathogenicity but also for informing the development of therapeutic strategies and vaccines aimed at combating current and future coronavirus outbreaks.

The article is organized as follows. In Sections 2 and 3, we describe the setup for multicanonical simulations used to investigate the folding dynamics of SARS-CoV RBMs and the microcanonical analysis framework employed to reveal structural phase transitions in folding, often missed by canonical ensembles. In Section 4, we present our results, focusing on the structural phase transitions identified through density of states calculations and the application of the Poisson-Boltzmann equation to model RBM solubility. We also discuss the broader implications of these findings, particularly how RBM mutations influence viral evolution, pathogenicity, and potential therapeutic strategies. Finally, in Section 5, we conclude by summarizing our key insights and proposing directions for future research.

We perform massive parallel Monte Carlo simulations [49] within a computational setup previously employed on modeling Amylin, Insulin, and Amyloid-$\beta$ interactions [43, 48, 50]. Here, the focus is on simulating the receptor-binding motifs (RBMs) of different SARS-CoVs to explore their folding dynamics and thermodynamic properties. The simulations employ the ECEPP/3 force field and the implicit solvation model SCH2, as implemented in the SMMP 3.0 package [47], which provides an effective balance between accuracy and computational efficiency, consistent with Molecular Dynamics predictions [43, 51, 52]. While the absence of explicit solvation is known to introduce phase transitions at significantly elevated temperatures, this is a minor drawback that must be accepted as a tradeoff for the high-statistics data achieved. This approach involves up to 36 million molecular sweeps across 300 multicanonical (MUCA) recursions, enabling a robust microcanonical analysis.

It shall be noted that RBMs of SARS-CoV-1 and SARS-CoV-2 has not been experimentally characterized in isolation, as available structural data pertain to the RBM embedded within the receptor-binding domain (RBD) of the Spike protein. The absence of experimental data on the isolated RBM poses a challenge, as its structural and dynamic properties may differ significantly from those observed in the native Spike context. This implies the need of computational modeling to generate initial structures suitable for simulation studies.

To begin the simulations, the RBMs for SARS-CoV-1, SARS-CoV-2 and its VOCs must first be encoded in the one-letter FASTA format, shown in (Table 1). Noteworthy, the $\beta/\gamma$ VOCs share the same RBM, whose mutations E484K and N501Y are prone to immunity evasiveness [19, 53]. Then, these sequences need to be converted into three-dimensional input structures (Figure 1), for instance using the I-TASSER (Iterative Threading ASSEmbly Refinement) software for homology modeling [54, 55, 56]. I-TASSER is a robust tool for protein structure and function prediction. It operates by initially threading the query sequence through a database of known protein structures to identify appropriate templates. The software then constructs full-length models via iterative fragment assembly simulations. This process not only predicts the 3D structure of the proteins but also offers potential insights into their biological functions, binding sites, and interactions.

The Multicanonical simulations are carried out in a rigid cubic box with 200Å sides, and the Boltzmann constant is set to $k_{B}=1.987\times 10^{-3}$ kcal/mol/K. Although initial molecular configurations are required to start the simulations, after thermalization, thermodynamic properties become independent of these initial conditions, as observed in ergodic Monte Carlo simulations. This allows for a comprehensive analysis of the folding dynamics and solubility of the RBMs across a wide temperature range, spanning from the homogeneous nucleation of ice crystals at 224.8 K to the critical point of water at 647 K.

Canonical Monte Carlo (MC) simulations are limited in their ability to explore the full energy landscape of complex systems due to their reliance on the Boltzmann weight [27, 30], which is defined at a fixed temperature. This constraint often restricts the utility of re-weighting techniques for extrapolating thermodynamic quantities across different temperatures [57]. In contrast, Multicanonical (MUCA) simulations [24, 25, 26, 27] offer a more versatile approach by employing generalized weights that allow the system to sample a wide range of energy states uniformly. The MUCA weight is given by

where $\Omega(E)$ denotes the density of states (DoS), $S(E)$ is the microcanonical entropy, $\bar{\beta}(E)$ is the microcanonical inverse temperature, and $\bar{\alpha}(E)$ is a dimensionless free energy term.

The MUCA method comprises two key steps: (1) determining the appropriate multicanonical weights, usually through an iterative process that involves histogramming the internal energy, and (2) performing a Markov Chain Monte Carlo (MCMC) simulation using the weights from the first step, followed by reweighting to the Gibbs ensemble. In this study, we focus on microcanonical thermostatistics to analyze the folding and solubility of the RBM from different SARS-CoV variants. Thus, our objective is to concentrate on step 1, specifically on computing the multicanonical entropy

where $\left\{\beta_{k},\alpha_{k}\right\}$ are piecewise functions representing the MUCA parameters, which approximate the microcanonical entropy.

The calculation is implemented through an iterative process [25, 27] where initial MUCA weights are set to $\omega_{muca}^{(0)}=1$, and the system undergoes a conventional Metropolis simulation. Data collected during this simulation is used to construct histograms of the energy distribution, $H(E)$, which are then used to update the MUCA weights. The iterative refinement of these weights follows the relation

where $c$ is a normalization constant ensuring that $S_{muca}^{(n+1)}(E)$ correctly represents the microcanonical entropy. Traditionally, iterations aim for histogram flatness; here, we instead monitor for stabilization of $S_{muca}(E)\times E$ in the whole energetic range, as described in [48], which takes $\mathcal{O}\left(5M\right)$ sweeps.

After convergence, the microcanonical thermodynamics [28] of the RBMs folding is automatically analyzed using the PHAST package [58]. The temperature $T(E)$ as a function of energy is derived from the microcanonical inverse temperature

where $k_{B}$ is the Boltzmann constant. The microcanonical specific heat $c_{V}$ is then defined as

Finally, the Helmholtz free energy $F(E)$ can be obtained via a Legendre transform at a fixed temperature $T_{c}(E)$

We thoroughly investigate for potential first-order structural phase transitions involved in the folding of SARS-CoVs RBMs. In such cases, microcanonical caloric curves $\beta(E)^{*}\times E$ often exhibit backbends, or S-bends, which reveal metastable states that remain hidden in canonical ensemble simulations [28, 31, 34]. The pseudo-critical inverse temperatures $T_{c}(E)$ associated with these transitions are determined using the Maxwell construction across the energy range $\Delta\tilde{L}=E_{fold}-E_{unfold}$, where $\Delta\tilde{L}$ represents the latent heat of the folding transition, and $E_{fold}$ and $E_{unfold}$ correspond to the internal energies of the folded and unfolded states, respectively [31].

Furthermore, the solubility of RBMs in water, which is critical for their biological function, is analyzed using the Poisson-Boltzmann equation [40, 48, 52]. The solvation free energy $\Delta G_{solv}$ is decomposed into nonpolar ($\Delta G_{np}$) and electrostatic ($\Delta G_{elec}$) contributions

The nonpolar contribution, associated with the formation energy of the molecular cavity, is often calculated using the solvent accessible surface area (SASA) model which is approximated as $\Delta G_{np}=\gamma\,\text{SASA}$, where $\gamma$ represents the surface tension, typically treated as an atom-independent parameter [59]. The electrostatic contribution is given by

where $q_{k}$ denotes the charge of atom $k$, and $\Delta\varphi_{vs}(r_{k})$ represents the electrostatic potential difference between vacuum and solvent environments at the position of atom $k$. Accurately solving the Poisson-Boltzmann equation (PBEQ) to obtain these potentials is computationally demanding but necessary for precise estimation of the RBMs solvation free energy. In this study, we use the PBEQ-Solver [40, 41, 42] for these calculations, with the dielectric constants set to $\varepsilon_{s}=78.5$ for the solvent and $\varepsilon_{p}=2.0$ for the protein, and the salt concentration set to 150 mM, as in Ref. [48].

Parallel simulations were performed to determine the MUCA parameters $\left\{\beta_{k},\alpha_{k}\right\}$, which were then utilized as input for the PHAST package [58]. These parameters enabled the derivation of the microcanonical thermostatistics governing the folding structural phase transitions of the RBMs in SARS-CoV-1, SARS-CoV-2 WT, and its early VOCs $\beta/\gamma$. Thermostatistical results are summarized in extensive physical units in (Table 2), they provide insights into the folding temperatures ($T_{fold}$), solvation energies ($\Delta G_{solv.}$), latent heats ($\Delta\tilde{L}$) and free energy barriers ($\Delta F$) that dictate the structural stability and phase transition behavior of these motifs. The caloric curves and intensive energetic barriers ($\Delta f=\Delta F/$residues) of the Helmholtz free energy at the (pseudo) critical inverse temperature ($\beta_{c}^{*}$) are depicted in (Figure 2) as functions of the internal energy per residue ($\varepsilon=E/\text{residues}$). These plots highlight the existence of regions with canonical metastabilities during the RBMs folding transitions, where using Maxwell’s constructions we identified the critical $\beta_{c}^{*}$ and the associated intensive latent heats ($\Delta L=\Delta\tilde{L}/$residues) for the SARS-CoV-1 and SARS-CoV-2 WT RBMs. In contrast, the $\beta/\gamma$ variants exhibit a second-order phase transition, marked by the absence of latent heat and Helmholtz free energy barriers, indicative of a smooth, continuous folding process similar to that observed in intrinsically disordered proteins such as Amylin and Amyloid-$\beta$ [43, 48, 50]. The critical $\beta_{c}^{*}\left(\varepsilon_{c}\right)$ of folding for $\beta/\gamma$ variant is detected by locating the peak of its $c_{v}\left(\varepsilon\right)\times\varepsilon$ (Eq. 5) — not shown. Furthermore, (Fig. 3) presents the isoelectric surfaces of representative RBM configurations for the viruses, sampled at physiological temperatures ($\sim 310$ K), which were used as input for computing the solvation energies via the PBEQ solver.

The correlation between thermostatistical outcomes and thermo-structural features during RBM folding transitions is crucial, especially in the absence of experimentally resolved RBM structures. Our MUCA simulations, quantifying residues in $\alpha$-helical and $\beta$-sheet conformations, illustrate how simulations might provide valuable guidance for future experimental design. Figure 4 shows the isocontours of internal energy per residue $\left(\varepsilon\right)$ as a function of residues in these configurations. Residue classification follows the SMMP framework [47], using dihedral angles $\left(\phi,\psi\right)$ with ranges of $\left[-70^{\circ}\pm 30^{\circ},-37^{\circ}\pm 30^{\circ}\right]$ for helices and $\left[-150^{\circ}\pm 30^{\circ},-150^{\circ}\pm 30^{\circ}\right]$ for sheets. This figure highlights two notable findings: SARS-CoV-2 is more $\alpha-$helix-prone than SARS-CoV-1, and the E484K and N501Y substitutions not only shift the folding phase transition to second order but also significantly enhance the $\alpha-$helical content in the $\beta/\gamma$ variants’ ground state compared to the WT. This trend, reminiscent of amyloidogenic Amylin isoforms studied using similar methods [58, 51, 52], can be related to the increased pathogenicity of these variants, as the Spike protein domains are known to induce protein aggregation under suitable biological conditions [44, 60, 61].

Next, we delve into the specific thermostatistical and solvation properties of each RBM, providing a detailed analysis and discussion.

This study provides a detailed microcanonical analysis of the RBMs in SARS-CoV-1, SARS-CoV-2 WT, and the $\beta/\gamma$ VOCs, offering critical insights into how thermodynamic properties influence viral evolution and pathogenicity. By utilizing the ECEPP/3 force field with implicit solvation, the analysis effectively balances computational efficiency with molecular accuracy, enabling robust data collection despite trade-offs such as elevated structural transition temperatures. The results highlight distinct thermodynamic and structural differences across these RBMs: SARS-CoV-1’s RBM, with a strong first-order phase transition and high solvation energy, exhibits structural rigidity that limits mutational flexibility and impedes immune evasion, reducing transmissibility. In contrast, the SARS-CoV-2 WT RBM shows greater thermodynamic flexibility, marked by a weaker first-order phase transition and reduced solvation energy, which promotes a broader mutational spectrum, facilitating the emergence of variants with enhanced immune evasion and transmissibility under selective pressures. The $\beta/\gamma$ VOCs, including the E484K and N501Y mutations, illustrate these evolutionary dynamics, increasing RBM solvation energy and shifting folding transitions to second-order, similar to intrinsically disordered proteins. This structural adaptation supports improved immune evasion and enhanced ACE2 receptor binding. Furthermore, the analysis reveals that SARS-CoV-2 RBMs exhibit a stronger propensity for $\alpha$-helix formation compared to SARS-CoV-1, with the $\beta/\gamma$ variants showing significant increases in $\alpha$-helical content in their ground states. This trend, reminiscent of amyloidogenic proteins like Amylin, may be linked to the increased pathogenicity of these variants, as the Spike protein domains are known to promote aggregation under certain biological conditions. These findings underscore the importance of ongoing surveillance and targeted interventions to mitigate the risks posed by emerging variants, as these thermodynamic properties suggest continued potential for viral evolution under selective pressures.

The Brazilian National Laboratory for Scientific Computing (LNCC) is acknowledged by supercomputing time granted at the Santos Dumont facility under project PHAST2. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. This article is a tribute to the memory of my esteemed colleague and friend, Dr. Marcelo Fernandes, who passed away prematurely due to COVID-19.