The inactivation of SARS-CoV-2 has become a major objective as the spread of the beta coronavirus (β-CoV) has triggered a global health crisis with high social (more than 37 million infected people and more than one million deaths worldwide as of September 29 th ) and economic (global GDP contraction up to 5.2 % for the year 2020) impacts (World Health Organization; World Bank). Virus inactivation can be achieved using different strategies based on chemical, biological and physical treatments. Biocide chemical agents and surfactants are frequently used to disinfect inanimate surfaces (Kampf et al. 2020a; Smith et al. 2020 ). Furthermore, virus in culture media are deactivated by chemicals as TRIzol ® , Formalin (formaldehyde) and β-propiolactone (Jureka et al. 2020) . Biological treatments, in the format of vaccines, are the most effective to fight the virus in living organisms, even though their successful development is sometimes limited by a combination of economic factors, regulatory environment and the empirical nature of modern vaccine discovery (Loomis and Johnson, 2015; Rappuoli et al. 2014; Tannock et al. 2020 ). Among physical treatments, cold plasma (Filipić et al., 2020) , far-UVC light (Buonanno et al., 2020) and membrane filtration (Shirasaki et al., 2017) have been described to eliminate virus in surfaces, air and water, respectively, even though thermal inactivation is the most extensively used treatment when possible (Abraham et al., 2020; Cimolai, 2020; Jureka et al. 2020; Kampf et al., 2020b; Pastorino et al. 2020; Yap et al., 2020) .
Focusing on the inactivation by temperature, it has been reported that SARS-CoV-2 was highly stable at low temperatures (i.e. 4 °C up to 14 days; Chin et al., 2020) but it is sensitive to heat. Indeed, several minimal inactivation temperatures, which are comprised between 56 ºC (45 min for total inactivation) and 100 ºC (< 5 min for total inactivation; Jureka et al. 2020), have been reported (Abraham et al., 2020; Cimolai, 4 2020; Jureka et al. 2020; Kampf et al., 2020b; Pastorino et al. 2020; Yap et al., 2020) . On the other hand, the infectivity of SARS-CoV-2 in solution was strongly reduced (up to 100-fold), even though the virus remained similarly infective in surfaces at 4 and 30 ºC (Kratzel et al., 2020) .
Temperature is known to promote changes in the molecular structure of biomacromolecules (i.e. nucleic acids, proteins, lipids) until affecting their functionality. As part of natural evolution, such alterations have been used by microorganisms as an adaptation mechanism to environmental changes, for instance by variations in thermal labile hydrogen bonds that connect the strands of nucleic acids and proteins (Sengupta and Garrity, 2013) . In the particular case of proteins, temperature is known to induce variations of the secondary and tertiary structure, causing structural alterations that modify their stability and their role in regulating cellular processes, signal transduction and intrinsic enzymatic properties. However, the response of proteins to changes in the temperature conditions can be very different. For example, some proteins have high thermal stability while others can unfold or even denature at moderate temperatures (Dong et al., 2018; Julió Plana et al., 2019; Lopes-Rodrigues et al., 2018) .
In this work we aim to unravel the effect of the temperature on the molecular structure of the SARS-CoV-2 spike glycoprotein, which is a homotrimer with three monomers with identical primary structure, named chain A, B and C. Coronaviruses use the spike to bind cellular receptors, triggering a cascade of events that leads to cell entry (Song et al., 2018; Zhou et al., 2020) . The spike protein is reported to bind the cellular receptor human angiotensin-converting enzyme 2 (ACE2) that mediates the fusion of the viral and cellular membrane and facilitates the virus introduction inside the cell (Ou et al., 2020) . Each spike monomer (180 kDa) contains 1273 amino acids (aa) and 5 consists of a signal peptide (aa [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] and two subunits, named S1 (aa and S2 (aa 686-1273), which are responsible for receptor binding and membrane fusion, respectively. The S1 subunit involves the N-terminal domain (NTD; aa and the receptor binding domain (RBD; 319-541). The latter consists of a core region with 5 pleated sheets (1, 2, 3, 4 and 7) organized in antiparallel model and the receptorbinding motif with two short -pleated sheets (β5 and β6), loops and alpha helices (α4 and α5; Shang et al., 2020) . A total of three cysteine residue pairs provide stability to the core and an additional cysteine residue pair helps in connecting the distal end of the receptor-binding motif . The S2 subunit has a key role in the membrane fusion (Kirchdoerfer et al., 2016) . Thus, after the initial interaction of the RBD of the S1 subunit and the peptidase domain of ACE2, the fusion with the host cell membrane continues through the interaction between the heptad repeat 1 domain (HR1; aa 912-984), the central helix (CH; aa 985-1035), the connector domain (CD; aa 1076-1141), the heptad repeat 2 domain (HR2; aa 1163-1213) of S2 to form a helix bundle fusion core (Jan Bosch et al., 2004) , which is extra stabilized by a short residue sequence (named fusion peptide of S2 or FP; aa 788-806; Xia et al., 2020) .
Furthermore, another important structural aspect of the SARS-CoV-2 spike is the identification of two states for the hinge-like conformation of the RBD of S1 from chain B (Walls et al., 2020; Wrapp et al. 2020 ). These are: 1) the closed down state ( Figure   1a ) in which the RBD of S1 covers the apical region of S2 near the C-terminus of HR1 (i.e. the receptor binding for interacting with ACE2 is buried); and 2) the open up ( Figure 1b ) state in which the RBD is dissociated from the central axis of S2 and the NTD of S1 (i.e. the receptor binding motifs are exposed). Moreover, this in silico methodology has been used to trace the atomic level contacts and interactions, which cannot be captured experimentally, and unveil a molecular mechanism for the virus inactivation. Such mechanism is intended to contribute to the development of new inactivation strategies by means of physical or chemical treatments specifically designed to disable such key interactions identified through the functional sites of the virus.
The atomic coordinates of the homotrimeric spike glycoprotein of SARS-CoV-2 in the closed down and open up conformational states were taken from the Protein Data Bank (PDB; entry 6vxx and 6vyb, respectively, from cryo-electron microscopy structures of the SARSCoV-2 spike ectodomain trimer; Walls et al., 2020) . Missing protein segments were built using the UCSF Chimera program (Pettersen et al., 2004) , which used the Modeller algorithm (Webb and Sali, 2016) to fill with the FASTA amino acids sequence while maintaining the crystallographic structure fixed. The different models generated for each conformational state were assessed using Z-DOPE (Discrete Optimized Protein Energy), a normalized atomic distance dependent statistical potential based on known protein structures. For each conformational state, the model 8 scored with the lowest Z-DOPE, which was lower than -1 for both cases, was selected for MD simulations. After adding the hydrogen atoms and forming the disulphide bonds, the structures were solvated, thermalized and equilibrated using the protocol described in the Methods section (Supplementary Information). Finally, 150 ns long MD production runs in an NVT ensemble were conducted at 298, 310, 324, 338, 358 and 373 K using the Amber18 program (Case et al., 2018) . In order to ensure the repeatability of the observed tendencies, replica simulations were performed at defined temperatures by changing the initial velocities. Accordingly, a total of 3.6 µs (150 ns × 2 conformational states × 6 temperatures × 2 replicas) were simulated for such a large system, which involved more than half million atoms.
The influence of the temperature on the closed down and open up states is quantitatively analyzed in Figure 1c -d, which compares the temporal evolution of two different parameters calculated using data derived from simulations at 298 K and 373 K.
These parameters, which have been defined to clearly differentiate between the two conformational states as well as to be able to identify the presence of intermediate Although detailed energetic analyses of the closed down and open up conformations are out of the scope of this work, we observed that the former state is stabilized with respect to the latter one, independently of the temperature. This stability order is in agreement with experimental observations (Pallesen et al. 2017; Walls et al., 2020; Yuan et al., 2017) . Even though the two conformational states are maintained throughout all the trajectories, temperature affects the spike structure, causing changes that, after a certain threshold, cause the inactivation of the virus (Abraham et al., 2020; Cimolai, 2020; Jureka et al. 2020; Kampf et al., 2020b; Pastorino et al. 2020; Yap et al., 2020) . These changes are illustrated in Figure 1e -f, which show the closed down and open up conformations obtained at the end of the MD trajectories at 373 K. As can be seen, the structures present significant differences to naked eye with respect to those shown in Figure 1a -b, which were used as a starting point. These differences affect not only the secondary structure, but also the quaternary structure, which describes the interchain assembly between the three monomers. Thus, visual
analysis of the open up conformation suggests that one of the monomers undergoes densification in certain regions due to structural deformations occurring at 373 K (blue monomer in Figure 1f ).
The influence of the temperature on the interchain assembly was analyzed by exploring the temporal evolution of the interchain distances (i.e. distance between two monomers), which was calculated as the distance between mass centres. Figure 2a In order to ascertain if the effect on the assembly of the chains is associated to a change in the global shape of the individual monomers and of the whole homotrimeric construct, the dynamics of the radius of gyrations (R g ) was followed. The R g for the individual chains does not exhibit abrupt changes when the temperature increases from 14 In order to look in detail the effect of the temperature on the different domains for the two states, Figure 4 compares the RMSF calculated at 298 and 373 K considering all atoms of each residue of chains A, B and C. As is shown, the position of the following domains is indicated for each subunit of each chain: 1) the NTD and the RBD from the S1 subunit; and 2) the FD, the HR1, the CH and the CD from the S2 subunit. For chain A, the S1 subunit is much more affected by the temperature than S2, reflecting that sheet rich domains, such as the NTD and the RBD, are more flexible than the helical rich domains located at S2. The flexibility of S1 is higher for the open up than for the closed down, as is evidenced by the fact that the RMSFs values reached at 373 K are, in general, higher for former state than for the latter one. Instead, the rigidity of S2 is similar for the two states with exception of HR1 domain, which shows large fluctuations at the central region for the closed down state. Inspection of the RMSFs achieved at intermediate temperatures indicates that the fluctuations reached by the S1 subunit of chain A increases progressively with heating. This is illustrated in Figure S1 , which displays the RMSF calculated at 324 and 358 K The behavior shown by chain B is similar to that of chain A (Figures 4 and S1 ), the S1 subunit being much more influenced by the increase in temperature than S2.
However, it should be noted that the fluctuations observed at 373 K in the NTD and RBD domains are much higher for the B chain than for the A one, indicating that in this case the flexibility of the S1 subunit is greater. This increase in the flexibility of S1 is particularly striking for the open state that reaches RMSF values close to 12 Å (NTD) or even higher (RBD). The effect of temperature in the open state is even more pronounced for the C chain, which shows a marked destabilization of the NTD at 373 K. It is worth noting that the thermal disruption of the NTD for the open up 15 conformation also occurs at 358 K ( Figure S1 ), while the RSMF profiles of the two states at 324 K are similar to that calculated at 298 K. In order to assess how the temperate affects the network of hydrogen bonds, a detailed analysis on this specific interaction was conducted. The geometric parameters and the cut-off values used to define the D-H···A interaction (where A is an acceptor atom, D a donor heavy atom, and H a hydrogen atom) as hydrogen bond are both the DHA angle, which must be greater than 135º, and the D···A distance, that must be lower than 3.0 Å. Figure 
The impact of the restrictions on the protein conformation has been further analyzed by considering distortions in the backbone ,-rotamers. For this purpose, ,-Ramachandran plots have been depicted for selected snapshots at the beginning (t= 0 ns), middle (t= 75 ns) and end (t= 150 ns) of the simulations. The approximate location and shape of the most favorable low-energy regions in this coordinate system is known to depend on the chemical structure of the residue. In this work, we have applied the classification proposed by Richardson and coworkers (Lovell et al., 2003) , according to which ,-maps of the following residues should be considered separately: 1) glycine (Gly) that is the most flexible residue due to the lack of side chain; 2) proline (Pro) with a singular cyclic structure that prohibits the rotation about the N-C  bond and confines the  torsion angle at around -60º; 3) residues preceding Pro (pre-Pro), which exhibit a very distinctive ,-distribution due to the steric restrictions imposed by the neighbor The ,-maps of General residues for the two states, calculated at 298 and 373 K, are shown in Figure S7a (t= 0 ns) and 7a (t= 75 and 150 ns). It is worth noting that the main part of such residues is located inside favored and allowed regions (blue dots within the contoured zones), while a small part is out but surrounding the allowed regions (red dots over the border or close to the contoured zones). The latter residues, which are found in a similar number for the two states, are modestly strained with a small energetic penalty with respect to the residues inside the favored regions. In addition, a few residues appear at forbidden areas with pronounced steric clashes (red dots over white regions), indicating highly strained conformations. These few disfavored residues appear in the closed state even at the beginning of the simulation ( Figure S7a) , independently of the temperature, suggesting that they should be associated to the symmetric assembly of the three monomers rather to the heating process. Moreover, highly strained conformations are localized in chains B and C, as is proven in Figure S8 that compares the ,-maps for the three monomers at the end of the 298 and 373 K trajectories.
This effect is less pronounced for flexible Gly residues, as is shown in the ,-maps displayed in Figures S7b and 7b . Gly lacks of a side chain, making its , values substantially less restricted than those related to other amino acids. Thus, the number of Gly residues localized in highly disfavored conformational regions (red dot over white regions) is very small, regardless of temperature. In this case, the increase in temperature from 298 to 373 K only causes a greater spreading of the blue dots inside the contoured regions, which is associated to the favored and allowed regions. This increment in the flexibility of Gly residue can be associated to structural deformations without significant energy penalties. A similar effect is observed for Pro and, especially, pre-Pro residues, as is shown in Figures S7c-d and S9 . Although temperature induces the apparition of a few strained residues for both the closed down and open up states, this straining effect is much less pronounced than for General residues. Besides, the apparition of strained pre-Pro residues is practically undetectable through all trajectories. It is worth noting that, in opposition to Gly, that is very flexible and able to accommodate in many different conformations, the behavior of Pro and pre-Pro has been attributed to the geometric restrictions imposed by the pyrrolidine ring. 
The stabilization of the closed down conformation with respect to the open up, which is in agreement with experimental observations (Pallesen et al. 2017; Walls et al., 2020; 23 Yuan et al., 2017) , has been attributed to the attractive inter-monomer interactions since the internal energy of each monomer is similar for the two states. patterns. More specifically, the higher fluctuations observed for the former state are more localized and less frequent than for the latter state, which exhibits important deviations at larger regions of the protein motifs. This feature is clearly illustrated by RMSF 298-373 values (Figure 3e-f ).
In general, the RMSFs profiles shown in Figure 4 are fully consistent with cryoelectron microscopy observations of SARS-CoV-2 spike, which revealed considerable flexibility and dynamics in the S1 subunit (Kirchdoerfer et al., 2016; Hoffmann, et al., 2020) . In addition of the expected flexibility of RBD, which has been extensively investigated in this work by considering the closed down and open up states separately, the NTD flexibility suggested by cryo-electron microscopy structures is confirmed by MD simulations. Another region from S1 that deserves consideration is the one located between the RBD and the S2 subunit, which is usually subdivided in two structurally conserved subdomains identified as SD1 and SD2. These subdomains, especially SD1, act as a hinge point for RBD closed down-to-open up transitions and, therefore, are strongly affected by the breathing movements between the two states, which explain the relatively large fluctuations observed at all studied temperatures.
On the other hand, the preservation of the number of hydrogen bonds with increasing temperature suggests that the glycoprotein spike of SARS-CoV-2 presents some kind of thermal stability, which at a first glance would not justify the inactivation of the virus by temperature. However, further understanding of the latter observation has been achieved by examining the changes in the topology of the hydrogen bonding network induced by the temperature. Thus, a large number of newly formed hydrogen bonds is observed at 373 K. Conversely, native hydrogen bonds were mostly preserved at 298 K for both studied conformational states through the dynamics runs. This difference explains the inactivation of the SARS-CoV-2 when rising the temperature. Although the largest change in the topology is achieved at 373 K, the observation of changes at lower temperatures is probably hindered by the limited length of the trajectories. More specifically, the total inactivation of the virus has been reported at 56 ºC after 45 min heating (Jureka et al., 2020) , which suggests a very slow kinetics for the conformational and topological changes related to hydrogen bonding reorganization.
Comparison of the topology maps at 298 and 373 K (Figure 5a-b) suggests that, in general, the newly formed hydrogen bonds involve residues that are relatively close to those involved in the native interactions. This feature is consistent with the fact that energy delivered through heating induces changes in the secondary and, even, tertiary structure of the spike glycoprotein but not in its global shape or quaternary structure, as discussed above (Figures 1 and 2) . Considering the moderate strength of hydrogen bonds (i.e. around -5 kcal/mol), the roughly retention of the shape at 373 K should be related to stronger noncovalent interactions. In fact, the average amount of salt bridges increases 11% (closed down) / 18% (open up) when the temperature rises from 298 to 373 K, supporting that salt bridges are responsible for the global shape retention observed at 373 K.
Overall, analyses of the temperature influence on both hydrogen bonds (weak directional interactions) and salt bridges (strong non-directional interactions) allow us to propose a molecular mechanism for the thermal inactivation of the virus. Such mechanism consists on differentiating the biological stability from the structural stability. The preservation of a large number of salt bridges and the formation of new ones upon increasing temperature are responsible of the structural stability in terms of protein global shape. Conversely, the deterioration of the hydrogen bonding network appears as the main reason for the functional inactivation of the virus upon heating. This mechanism is summarized in Scheme 1, which depicts the drastic conformational changes experienced by the NTD and RBD of chain B when temperature increases from 298 K to 373 K, whereas the overall shape of the spike is roughly preserved. The ,-maps displayed in Figures 7 and S7-S9 are fully consistent with the mechanism proposed in Scheme 1. Thus, the General residues, which are able to form hydrogen bonds and salt bridges, experience not only the largest conformational changes but also, in some cases, adopt strained conformations, as suggested by the topology maps calculated for those interactions. Such findings, open new avenues to develop strategies to inactivate the virus, by targeting specific areas of the homotrimer in order to misbalance, remodel or cleave the hydrogen bonds that make the virus functional (i.e. by using chaotropic agents or surfactants that make feasible the disruption of hydrogen and salt bridges). Furthermore, it allows the development of active devices able to deliver energy that can cleave and reorganize hydrogen bonds, for instance, through ultrasounds, irradiation or by using nanoparticles as nanosources of heat able to induce a local increment of temperature by using the local surface resonant plasmon effect, a phenomenon that occurs when light interacts with matter at specific wavelengths. The effectiveness of such local thermal treatment could be increased by controlling the intensity of the light source, the irradiation time lapse and the functionalization of nanoparticle designed to maximize the interaction with the virus areas where the higher number of hydrogen bonds is exposed.
In summary, our results suggest that temperature induces conformational changes on the S1 subunit of the spike glycoprotein of SARS-CoV-2 that remodel the internal hydrogen bonding structure. However, the impact of temperature, even at 373 K, is not strong enough to induce a significant modification in the global shape of the spike.
Those conformational changes, are much more pronounced in the state where the binding domain is accessible and the virus is infective (open up state) than when the 28 binding domain is retracted (closed down state). This effect has been associated to a drastic modification in the hydrogen bonding topology, which particularly affects the recognition functionality of the receptor binding domain. Such network reorganization is triggered by the energy delivered through heat that allows the massive cleavage and remodeling of a significant amount of hydrogens bonds. Conversely, the salt bridges topology remains much less altered, allowing to maintain the main structure that defines the shape of the spike. The deep knowledge about such inactivation mechanism facilitates the development of new strategies intended to inactivate the virus through the destruction or modification of such hydrogen bonds by chaotropic agents and surfactants or through physical treatments able to selectively target such labile hydrogen bond structures. The missing residues (i.e. 44-55, 88, 89, 118-139, 147-159, 217-236, 417-429, 445-462, 476, 595-614, 651-653, 802-828, 1162-1175, 1236-1258, 1268-1280, 1292-1294, 1307-1309, 1338-1357, 1550-1556, 1562-1585, 1610-1616, 1716-1735, 1772-1783, 1923-1948, 2283-2296, 2360-2380, 2389-2401, 2459-2479, 2661-2263, 2671-2677, 2687-2706, 2837-2856, 2893-2905 and 3044-3071) were incorporated using the Modeler algorithm (Webb and Sali, 2016) implemented in the UCSF Chimera program (Pettersen et al., 2004) randomly Na + and Clions until reaching a 0.15 M NaCl concentration to reproduce physiological conditions. Then, the two models were processes with the LEaP program (Case et al., 2005) to add hydrogen atoms to the protein and to generate Amber topology files and coordinates files. Accordingly, the models used to represent the closed down and the open up states of the homotrimeric spike protein presented 516706 and 558601 explicit atoms, respectively.
All simulations were performed using the AMBER 18 simulation suite (Case et al., 2018) . Protein atoms were modeled using the Amber ff14SB force field (Maier et al., 2015) , the glycan atoms included in the cryo-EM coordinates were modeled using the Glycam06 force field (Kirschner et al., 2008) , and water atoms were modeled using the TIP3P force field (Jorgensen et al., 1983) . Equilibration calculations were started by relaxing the protein regions filled with the UCSF Chimera program (Pettersen et al., 2004) , which was achieved by applying the
BroydenFletcher-Goldfarb-Shanno quasi-Newton algorithm methodology to the new added residues meanwhile the rest of the system was kept frozen. Next, the whole system was submitted to 2500 steps of full conjugate gradient minimization to relax conformational and structural tensions.
The Langevin dynamics method (Izaguirre et al., 2001) was used to heat the system and to rapidly equilibrate its pressure and temperature. The relaxation times used for the coupling were 5 ps for both temperature and pressure. The temperature was increased from 0 K to 298 K using 60 ps simulation in the NVT ensemble, using an integration time step of 1 fs and keeping the pressure at 1.034 bar. Then, 1 ns in the NPT ensemble were conducted at 298 K to relax the structure and the density (integration step: 2 fs; pressure:1.034 bar). The last snapshot of this relaxation was used not only as starting point of the NPT production trajectory at 298 K but also as starting point of 0.25 ns NPT-MD simulation to bring the temperature to 310 K (integration step: 2 fs; pressure:1.034 bar). The last snapshot was used as starting point for the NPT production trajectory at such temperature and as starting point for the NPT-MD used to bring the temperature to 324 K. The same process was used to generate starting points for the production simulations 338, 358 and 373 K.
Production NVT-MD trajectories at 298, 310, 324, 338, 358 and 373 K were 150 ns each. Replica simulations were performed at all the indicated temperatures by changing the initial velocities. Accordingly, a total of 3.6 µs (150 ns × 2 conformational states × 6 temperatures × 2 replicas) were simulated for the studied system.
Non-bonding pairs list was updated every 10 steps. Periodic boundary conditions were applied using the nearest image convention and the atom pair cut-off distance used to compute the van der Waals interactions was set at 10.0 Å. In order to avoid discontinuities in the potential energy function, non-bonding energy terms were forced to slowly converge to zero, by applying a smoothing factor from a distance of 12.0 Å.
Beyond cut-off distance, electrostatic interactions were calculated by using Particle
Mesh of Ewald (Toukmaji et al., 2000) . Figure S1 . Blue dots correspond to residues with allowed conformations (regions enclosed by the contours, see Figure S6 ), while red dots represent residues with strained (close to regions enclosed by the contours) or very strained (red dots over white regions)
conformations. 