Heat-induced transitions of an empty minute virus of mice capsid in explicit water: all-atom MD simulation

Abstract The capsid-like structure of the virus-based protein nanoparticles (NPs) can serve as bionanomaterials, with applications in biomedicines and nanotechnology. Release of packaged material from these nanocontainers is associated with subtle conformational changes of the NP structure, which in vitro, is readily accomplished by heating. Characterizing the structural changes as a function of temperature may provide fresh insights into nanomaterial/antiviral strategies. Here, we have calculated heat induced changes in the properties of an empty minute virus of mice particle using large-scale ≈ 3.0 × 106 all-atom molecular dynamics simulations. We focus on two heat induced structural changes of the NP, namely, dynamical transition (DT) and breathing transition (BT), both characterized by sudden and sharp change of measured parameters at temperatures, TDT and TBT, respectively. While DT is assessed by mean-square fluctuation of hydrogen atoms of the NP, BT is monitored through internal volume and permeation rate of water molecules through the NP. Both the transitions, resulting primarily from collective atomistic motion, are found to occur at temperatures widely separated from one another (TBT>TDT). The breathing motions, responsible for the translocation events of the packaged materials through the NP to kick off, are further probed by computing atomic resolution stresses from NVE simulations. Distribution of equilibrium atomistic stresses on the NP reveals a largely asymmetric nature and suggests structural breathing may actually represent large dynamic changes in the hotspot regions, far from the NP pores, which is in remarkable resemblance with recently conducted hydrogen-deuterium exchange coupled to mass spectrometry experiment. Communicated by Ramaswamy H. Sarma


Introduction
About half of the known virus families, the protein shell that encloses the viral genetic material, is a 'spherical' or icosahedral capsid (Crick & Watson, 1956), often ranging in diameter from 30 to 100 nm. These capsids are self-assembled systems composed of multiple copies of individual proteins. Assembly, stability, and disassembly of virus particles is a marvel of nature that involves massive coordination of protein-protein and protein-nucleic acid interactions. These interaction forces are although much weaker than covalent bonding, the capsid can withstand tens of atmospheric internal pressure (Evilevitch et al., 2003;Smith et al., 2001) before they could be ripped open. Certain Capsids (for example, cowpea chlorotic mottle virus, both empty and containing the RNA genome) can even recover from deformations as large as 30% before rupturing . In the minute virus of mice (MVM), the internal stress provides the initial impetus for release of the confined genetic material, upon its binding to a receptor in the outer membrane of the host cell (Carrasco et al., 2008). Empty capsids without the genetic material are also found among natural viruses (Ansardi et al., 1996;Curry et al., 1997;Roizman et al., 1958). The remarkable stability of the protein-based capsid-like nanoparticles (NPs) that can withstand high hydrostatic pressure and hazardous chemical environment are also ideal candidates to serve as nanomaterial for a range of novel applications such as nanocontainers in medicine and materials science (Douglas & Young, 1999;Kasuya et al., 2005;Koutsky et al., 2002;Lowy & Schiller, 1999;Patterson et al., 2012Patterson et al., , 2014Porta et al., 2013;Rome & Kickhoefer, 2013;Servid et al., 2013;Singh et al., 2006). Viral capsids are also guiding the design of engineered protein NPs to enhance mechanical strength against breakage (Medrano et al., 2019).
The release of the packaged material from these nanoscale machineries, such as in MVM, are associated with a subtle conformational change of the structure (Castellanos et al., 2012;van de Waterbeemd et al., 2017). In vitro, this change can be readily accomplished by heating, which helps in overcoming the conformational energy barrier and subsequent capsid disassembly, for example, in foot-and-mouth disease virus (Rinc on et al., 2014), and bacterial viruses (Qiu, 2012;V€ or€ os et al., 2018), Detection of the heat induced structural mobility and integrity of virus (Tresse et al., 2017) and virus-like NPs (VNPs) can help designing improved vaccines against thermal denaturation (Rinc on et al., 2014), or, antiviral agents (Lewis et al., 1998) Studying the successive structural changes as a function of temperature provides valuable insight into the principles of viral assembly (Hagan & Chandler, 2006), which in turn may help the inactivation procedure of viruses, vital for food industries and pharmaceutics (Hirneisen et al., 2010;Los et al., 2004;Muller-Merbach et al., 2005). In addition, the determination of whether, and to what extent, elastic properties vary across the VNP surface (Carrasco et al., 2008;Castellanos et al., 2012;Klug et al., 2006) as a function of temperature will be important in understanding the processes of virus disassembly. The aim of this work is precisely to investigate the heat induced changes in the microscopic properties of a VNP at the molecular level. Our central result is that the 'breathing' motion associated with the release of the packaged materials in the nanocontainers (viral infectivity), or, the disassembly of VNP lie in regions far from its pores. We also found that the mapping of the equilibrium atomistic stress on the VNP structure is a perfect mimic for the recently conducted hydrogen-deuterium exchange coupled to mass spectrometry (HDX-MS) experiment, aimed at delineating similar heatinduced structural changes (van de Waterbeemd et al., 2017).
A large number of experimental techniques have been applied in the solution phase to probe the capsid dynamics that includes atomic force microscopy (AFM)-based nanoindentation measurements (Hernando-Perez et al., 2014;Ivanovska et al., 2011;Klug et al., 2006;Kol et al., 2006;Medrano et al., 2019;Michel et al., 2006), radiolabelling of virions (Rinc on et al., 2014), differential scanning calorimetry measurements (Qiu, 2012), fluorescence measurements (Carreira et al., 2004;Poian et al., 2002), Raman spectroscopy (Verduin et al., 1984), HDX-MS (Lanman et al., 2004;Monroe et al., 2010;van de Waterbeemd et al., 2017;Wang & Smith, 2005), and limited proteolysis combined with peptide mass mapping by MALDI mass spectrometry (Bothner et al., 1999;Lewis et al., 1998;Speir et al., 2006). The last two techniques have proven to be very promising for identifying the dynamic regions of capsid proteins and also to measure relative differences in dynamic character of the peptide fragments. Molecular dynamics (MD) simulations (Hagan, 2014), on the other hand, are now routinely applied to study large biological systems, such as, virus capsids, at both coarsegrained (Grime et al., 2016;Nguyen et al., 2007;Ross et al., 2018) and all-atom levels (Andoh et al., 2014;Jiang et al., 2015;Larsson et al., 2012;Perilla et al., 2016;Perilla & Schulten, 2017;Zink & Grubmuller, 2009. Because MD simulations can provide information at the resolution of individual particles, the technique can be utilized to unearth the properties of nanocontainers in a much greater detail at the atomic level. Here we endeavor to apply MD techniques, at the atomistic level, to study a host of properties of an empty MVM particle as a function of temperature. MVM belongs to the family of parvovirus, the capsids of which are % 26 nm in diameter and contain 60 copies of viral protein (VP) 1 to 3 (VP1-VP3) in a T ¼ 1 icosahedral arrangement encapsulating single-stranded DNA (ssDNA). The parvovirus capsid surface is classified with protrusions at or surrounding the icosahedral 3-fold axes and depressions at and around the icosahedral 2-fold and 5-fold axes, respectively (Kontou et al., 2005). The N-terminal (Nt) domains of VP2 and VP1 are initially buried in the capsid, which are responsible for molecular signals that are required during various stages of the infection cycle, viz. packaging and release of the viral ssDNA. They are externalized through capsid pores located at the 5-fold symmetry axes. VP2 capsid proteins are also able to self-assemble into DNA-free empty VNP, which is structurally indistinguishable from natural MVM capsids. This VNP is an ideal model system to gain insight into the MVM capsid assembly/disassembly mechanisms. By construct, it facilitates specific structural investigation on the externalization of VP2 Nts, which is triggered in vivo by the confined DNA (Cotmore & Tattersa, 2014) and in vitro by heating at ! 319 K (Hernando et al., 2000;van de Waterbeemd et al., 2017).
As we demonstrate below that around this temperature, the capsid changes from a nonfunctional form to a flexible and functional state, such that, the breathing motions of VNPs become large enough to initiate protein translocation through capsid pore. We formally call this heat induced structural transition as breathing transition (BT) and denote the associated temperature as T BT . When combined with the atomistic mechanical stress data, our MD results can be found to be consistent with a dynamic or breathing capsid in which the different peptide fragments experience different tensile stress and specific capsid regions undergo massive structural changes, far from the capsid pore, located at the 5-fold symmetry axes. BT occurs primarily from the collective heavy-atom motions.
Like in a mature capsid, at sufficiently low temperatures, the biological activities of its isolated protein subunit also cease to exist. With increase in temperature they come into 'life', typified by a sudden and sharp transition to a nonlinear temperature dependence of certain measurable parameters, such as, atomistic mean square displacement (MSD) of the lightest H-atoms. In parallel with the glass transition temperature in a glassy substance, the onset temperature of anharmonicity of an isolated protein molecule is referred to as the dynamical transition (DT) temperature, T DT (Liu et al., 2017;Pathak & Bandyopadhyay, 2019;Yoon et al., 2014). DT is also expected to occur in a full capsid structure. Both DT and BT are the causal effect of the collective dynamics of the structure, which we quantify here by heat induced changes in the MSD of hydrogen atoms, internal volume, rate of water permeation through the semi permeable surface, and the elastic properties of the VNP. We will show below that for the VNP, T DT is much lower than T BT : while DT brings forth the flexibility in the VNP, BT brings it in a functional state, wherein heat induced structural rearrangements of the structure can be thought to be sufficient for the associated translocation events through the capsid pore to kick off.

Results and discussion
A large number of isoenergetic substates in the conformational energy landscape of a bio-macromolecule is generically similar to those in glass formers. Upon supercooling both can be trapped in metastable states with exceedingly large relaxation times. Consequent to heating, the activities of the macromolecule begin to emerge and beyond a certain onset temperature their atomistic motions show anharmonic temperature dependence. In line with glass transition in glassy material, it is tempting to associate the onset of anharmonicity with biomolecular glass transition. The decisive work by Stanley et al. (Kumar et al., 2006) has brought forward a paradigm shift in the scientific thought processes pertaining to anharmonic motion in a biomolecular system that has laid to present trait to identify the phenomenon as DT at the onset temperature, T DT (Liu et al., 2017). A nanocontainer, like a viral shell, is also expected to exhibit similar thermal behavior, although the corresponding temperature may not be enough to reach the threshold for its heat-induced structural breathing motions. Here we set out to investigate the two heat induced transitions, both caused by atomistic motions and show that like DT, the onset of BT at T BT can also be identified by sharp changes in certain physical parameters, and that T DT and T BT are widely separated.

Dynamical transition
The incidence of DT at $ 180-200 K, at which a hydrogenated protein chain acquires sufficient flexibility, so as to be detected by incoherent neutron scattering and MD studies, is well documented. However, little is known whether similar heat-induced nonlinear response can also be observed in protein assemblies. The participation of the heat bath through the solute-solvent coupling is irrefutable in the phenomenon (Gupta et al., 2016;Pathak & Bandyopadhyay, 2019). In a hydrogen bond network forming solvent (water), the collective dynamics of solute-solvent motions can also influence the internal motions of a VNP. The resulting temperature threshold for the sudden surge in atomistic fluctuations could be of considerable importance for materials applications of protein assemblies. Since H is the lightest element, its atomistic fluctuations in response to change in temperature of the surrounding heat bath will be first to be noticed. In view of this, in Figure 1 the temperature response of H atoms' MSD are plotted (cf. Equation (1)). In accord with the previous studies on protein systems (Gupta et al., 2016;Kumar et al., 2006;Liu et al., 2017;Pathak & Bandyopadhyay, 2019;Yoon et al., 2014), our simulations confirm the DT for a viral shell. We have also conducted separate simulations for an isolated protein subunit and five subunits (forming a capsomer) of the VNP and the results of their MSD H variations are also shown in the figure.
Clearly, for all the three cases, the MSD profiles changes from harmonic to anharmonic at T DT , evidenced by the changes in slope of their thermal response. The amplitude of MSD for all H atoms is always larger than that of the amide-H atoms, which is in line with the expectation that the polar amide-H atoms are also engaged in intramolecular hydrogen bonding. Note that although the amplitude of MSD variations for all H and amide-H atoms are different from each other, they change their harmonic to anharmonic behavior at the same temperature, T DT . From a 'frozen' cryogenic state, as the temperature increases, any molecular motions require the surrounding solvent molecules to accommodate the concomitant changes. This renders the solute protein's conformational free-energy landscape to be slaved to the surrounding solvent. For one single protein chain (or for that matter their five assembled units in a capsomer) there are easily accessible excitations available in the energy landscape. Whereas for the large assemblies of protein chains in a viral shell, such conformational transitions may not be so easily accessible. Extensively large intramolecular interactions in viral shell may render them to be trapped in deeper potential wells, such that, the solute-solvent coupled motion sets into action at elevated temperatures. Thus it is likely that the observed T DT for the shell is larger (by $ 30 K) than either of its subunits ( Figure 1). Note that to determine the DT, we used linear regression of the computed MSD values from 100 to 200 K and from 200 K to the final higher temperatures and extrapolated them. The intersection between the straight lines for the harmonic and anharmonic behavior is considered to be taken as T DT . If we included slight variations in the number of data points in the low or high temperature segments during the linear regression analysis, we obtained a worse correlation coefficient with a variation in intersection temperature by ±5 K. Similar linear regression analysis protocol is adopted while evaluating T BT (for example, Figure 2).
As discussed below, large structural changes of the VNP associated with the breathing motions of the capsid in which translocation of internal proteins to the exterior of VNP surface occur (Lewis et al., 1998;van de Waterbeemd et al., 2017) is different from DT. And the event takes place at a temperature (T BT ) much higher than T DT . Note that like the VNP at higher temperature, the isolated protein chain, one or five, will also undergo conformational changes (Yoon et al., 2014). However, the MSD of H-atoms (cf. Figure 1) in an equilibrium MD run is a net outcome of the resulting changes.

Breathing transition
The VNP stability balances on a fulcrum (Medrano et al., 2019), as they must be sufficiently stable to withstand the physiochemical insult due to the surrounding medium and yet also efficiently initiate the release of packaged cargo. Under certain triggers such as receptor binding, change in pH, or elevated temperature, the non-trivial structural changes associated with the cargo release become feasible. Here we present results for heat-induced major structural changes. First we determine the breathing transition temperature and then discuss the results on the breathing motions of the VNP.

VNP internal volume
In response to the temperature variations, the internal volume of the VNP must also vary. Utilizing Richards' rolling probe method the internal volume change of the VNP is calculated. In Figure 2 the results for % change in VNP capsid volume with respect to that at 273 K is presented as a function of increasing temperature. The results points to the shrinkage of capsid size between 273 K and 318 K. Conceivably, structural rearrangements (Ceres & Zlotnick, 2002;Kegel & van der Schoot, 2004) occur in the capsomeric proteins in this temperature range, which results into an enhanced interaction hence to volume decrement of the virus particle. This is also consistent with the observation that the capsid becomes more compact between 273 K and 318 K, as can be inferred from the measured radius of gyration (R gyr ) of its Ca atoms (see the inset in the Figure 2, R gyr at 273 K is larger than that at 318 K). However, beyond 318 K the measured R gyr values increases as the temperature increases further and the capsid structure losses its compactness, indicating the heat-induced structural changes of the capsid. Sudden sharp rise in volume change at a temperature, T BT is regarded as the onset of breathing motions of the capsid's peptide fragments. T BT is estimated to be 318 K Figure 1. Temperature dependence of MSDs of H atoms' fluctuation at 100 ps for all H atoms (open symbols) and amide H atoms (filled symbols) of (A) isolated one protein chain forming the VNP, (B) isolated five protein chains surrounding a VNP pore, and (C) the VNP shell in full. The bottom panel shows the respective molecular structures (from left to right: monomeric protein chain, 5 monomer units forming a capsomer, and 12 capsomer units forming the full VNP). Therein the residue numbers 153 to 164 are presented in gray colored VDW presentation that constitute 12 symmetrically placed pores on the VNP surface. The error bars in MSD profiles, estimated by computing 10 separate MSDs in various 100 ps time windows of the trajectories are smaller than the plotting symbols' thickness. The changes in the MSD profiles of the two H atom classes are seen to vary beyond T DT , and T DT of the shell is higher than its constituent subunits. Figure 2. Percentage change in VNP capsid volume with respect to its volume at 273 K as a function of temperature. Sudden sharp rise in volume change is identified as breathing transition occurring at T BT , as indicated. Inset: compactness of the capsid shell as measured by radius of gyration (R gyr ) of C a atoms of the pre-equilibrated VNP structures at three temperatures are shown. R gyr is calculated using the relation, where m i is the mass and r i is position of the ith C a atoms with respect to the center-of-mass (COM) of the shell. from the intersection of the linear fits to the high and low temperature data.

Water permeation
The semi permeable VNP allows water to be exchanged between its interior and exterior. For determination of the interior versus exterior regions of the VNP capsid see Computational Details section below. Permeation rates for water is calculated at various temperatures from the number of water molecules that were located on one side of the VNP surface at a reference time t and found on the opposite side after Dt time is elapsed (cf. Equation (2)). The slope of the linear fit to the data determines the permeation rate. To check the consistency, the time windows for calculations are varied. In Figure 3 we have presented the results for permeation rate versus temperature variation along with the consistency check of the calculation at a representative temperature (318 K). In line with the enhanced compactness of the shell between 273 K and 318 K temperature range (see discussion on VNP internal volume above), water permeation rate also found to decrease from 217 molecules/ns at 273 K to 180 molecules/ns at 318 K. That is, for an average 100 000 water molecules inside the capsid, it takes about 460 ns at 273 K and 555 ns at 318 K for them to be completely exchanged with water from exterior. However, beyond 318 K, in parallel with the capsid volume, the permeation rate can also be seen to increase. This is indicative of the loss of the VNP's structural integrity at elevated temperature. The intersection of the linear fits to the high and low temperature data can be thought to be the onset of BT at 318 K.
Having determined the BT temperature, below we set out to investigate the principle peptide fragments of the VNP structure responsible for BT and identify the key breathing motions. Our choice of protein fragments is largely inspired by a similar study on heat-induced structural transitions through HDX-MS experimentation (van de Waterbeemd et al., 2017).

Breathing motion
Self-assembled viral particles are exceptional for their accurate synchronization of a large number of capsid proteins, which renders maximum enclosed volume to hold on to the packaged material inside. For successful release of the inside material, the viral assembly must undergo sufficient structural changes. For example, during endocytosis, the protease gains temporary access to the Nts of capsid subunits, which are located at the interior of the shell (Kontou et al., 2005). The breathing transitions, which involve translocation of internal proteins to the exterior of the viral surface may actually represent large dynamic changes in the shell structure (Lewis et al., 1998;van de Waterbeemd et al., 2017). In above, to explore the thermal lability of the VNP, macroscopic theories have been utilized that neglect atomistic details of the capsid. And hence they remained mute about the most dynamic temperature sensitive hotspot region of the shell. Here we endeavor to identify the hotspot regions by connecting macroscopic stress to microscopical forces though atomistic mechanical stress-strain calculations (Figures 4-6).
The key element of the icosahedron structure is that it contains fivefold symmetry axes. The assembly of such subunits into closed surfaces is a resultant of curvature, introduced by the interactions between the pentagons, enclosing space with twelve pentagonal faces. The icosahedron, thus economizes the surface area-to-volume ratio. In this sense the viral shell can be considered to be structurally perfect. And therefore it is expected that the intrinsic and fundamental physical properties of the shell  would determine the thresholds for structural disassembly under stress at elevated temperatures. Determination of the corresponding transition temperature is also of paramount importance for preserving medicinal biochemicals, such as, attenuated vaccines. For temperatures below this, it greatly . Water permeation through the semi permeable VNP capsid. (A) water permeation rate as the temperature is varied, which exhibits the breathing transition at T BT , and (B) accumulated number of waters at 318 K crossing the surface of the shell at different time windows toward out (red) and toward in (blue), defined by the number of waters that were located in one side of the capsid at t À Dt and found in the other side at a later time, t: The slope of the linear fit to the data determines the permeation rate going out and in, which is found to be 178 ± 9 and 180 ± 12 water molecules per ns, respectively, over several time windows of trajectory, as shown. The similarity of the two rates is a signature of the convergence of simulations.
reduces molecular mobility of the packaged genetic material (Pathak & Bandyopadhyay, 2017) in addition to maintaining the shell's structural stability.
Atomistic force distribution analysis of the shell structure could provide valuable insights unraveling the molecular mechanism behind the BT phenomenon. To achieve this  goal, calculations of stress tensor at temperature-varying heat baths are performed. The mechanical properties of the system as measured by the global stress can also be measured at the atomic level as the equilibrium average of the corresponding virial, r ij (cf. Equation (3)) (Basinski et al., 1971;Fenley et al., 2014;Hatch & Debenedetti, 2012;Zimmerman et al., 2004). In Figure 4 we have presented the results of calculated stress levels on the VNP's peptide fragments by computing the virial stresses on individual atoms at various heat baths. In the Supporting Information (SI), Table S1, complete list of all peptide fragments and the corresponding temperatures at which they attain maximum stress level are provided. The results call for an explicit attention to the inhomogeneity of stress distribution on the shell that is inherent to its discrete and polyhedral nature. As is well known that conformational stress often gives rise to high plasticity in virus, for example, in Hepatitis B (Bottcher et al., 2006). In particular, temperature mapping of the stress distribution at which the fragments attain their highest stress values (Figure 4) reveals: dissimilar nature of the stress distribution between inner and outer surface of the capsid wall. For example, peptides forming the inner wall attains highest stress level at an early temperature (Figure 4(B), green colored fragments at 290 K, which constitute a large fraction of the mapped capsomer). Whereas peptides lining the outer surface show such propensity at temperatures higher than 290 K. As we discuss below, the motions behind the thermal fluctuation of the capsid can be inferred from the thermal Figure 6. Slopes to attain maximum stress from temperature T i to T j are mapped onto five subunits surrounding a capsid pore (at center). Panel A and B show outside and inside view, respectively. For visual clarity only one capsomer structure is shown. The peptides that experience larger slope (> 15) are colored purple while for smaller slope (< 15) they are colored green. The temperature intervals (T i -T j ) are indicated at the top of the panels. (C) Slope change to attain maximum stress from T i to T j and then to fall from T j to T k for each peptide segments are shown for the temperature intervals T i -T j -T k , as indicated in the panels. Color code used: red for maximum slope at 318 K, purple for maximum slope at 323 K, black for maximum slope at 328 K, and orange for maximum slope at 333 K. response of the stress levels of its constituent protein segments, which may eventually help identifying the hotspots most actively involved in initiating the conformational change across the protein assembly during BT.
Both tensile and compressive stress may lead to VNP disassembly. In panel (A) of Figure 5, total stress accumulated in a 1 ps trajectory window for three peptide fragments (for the choice of these particular fragments see below) are shown as a function of temperature, while in panel (C) an animated view of the distribution of tensile and compressive stress in these fragments are presented in the relevant temperature segments. In Figure 5(B), the result for time-averaged total stress of the full VNP versus the seven temperature points considered in Figure 4 is presented. Sudden sharp change in stress-temperature profile can again be identified as BT occurring at 320.5 K, which is intimately close to T BT value (318 K) as obtained earlier through other macroscopic quantities (shell volume, water permeation rate). This resemblance reaffirms the close connection of macroscopic observables to microscopically atomistic stress (Zimmerman et al., 2004).
The externalization of N-termini from the VNP is orchestrate by translocation of peptides through capsid pores, located in the center of its 12-pentamer structures. HDX-MS experimental results had previously demonstrated that the peptide translocation event during BT is actually initiated by subtle conformational change of peptide regions away from the pores (van de Waterbeemd et al., 2017). Enlightened with this observation, we further set out to identify the hotspot regions (primarily responsible for the global conformational rearrangement of the shell structure) and seek comparison of results on the peptide fragments calculated on the basis of their atomistic stress with those inferred from HDX-MS experiment. In this regard we have measured the stress levels of each peptide fragments at temperatures 290, 310, 318, 323, 328, 333, and 400 K (see SI, Figures S1-S3 for complete list) and calculated the slope of the stress-temperature profiles between two temperatures (SI, Table S1 full list is provided).
Following the results presented in Figure 4 for the attainment of highest stress levels, we further classify the peptide fragments into four thermal-response categories: those in temperature ranges 310-318 K, 318-323 K, 323-328 K, and 328-333 K where their stress levels attain the maximum values and fall thereafter at a subsequent higher temperature. In analogy with deuterium uptake plots in HDX-MS experiment (van de Waterbeemd et al., 2017), a higher slope to attain the maximum stress level from temperature, T i to T j can be equated with higher dynamic mobility of a given peptide fragment. Similarly, a higher slope (negative) to fall from the maximum stress level later at a higher temperature, T k would imply dynamic mobility of the peptide fragment decreases abruptly. In Figure 6 we have summarized the results. The heterogeneous thermal-response of the peptide fragments are quite evident. Slopes of the stress-temperature plots (cf. Figure 5(A)) of each of the peptide fragments (see bottom panel, Figure 6) classified, for example, in the 318-323-328 K temperature segment is the highest while heating from 318 K to 323 K and their stress levels fall thereafter as temperature increases further to 328 K. This would correspond to transition temperature, 320.5 K (see for example, Figure 5(A) for the peptide fragments 178-195 and 248-255). Incidentally this temperature is also recorded in the total stress-temperature profile of the full shell (see Figure 5, panel B), which has been previously identified as T BT . Whereas peptides in the 310-318-323 K, 323-328-333 K, and 328-333-400 K temperature segments [ Figure 6(C)] either attain the maximum stress levels before or after the said transition temperature (i.e. they are either flexible enough at low temperature or attains sufficient tensile flexibility after the BT has been initiated). For example, consideration of fragment 196-215 (in the 323-328-333 K segment) would result in to higher T BT (325.5 K) ( Figure 5, panel A). It is important to emphasize here that given the approximate nature of force field calculations, the discovered transition temperature values for the three peptide fragments are identical with the HDX-MS experimental results (van de Waterbeemd et al., 2017). Although this could as well be a fortuitous outcome, the remarkable resemblance points to the close correlation between deuterium uptake and atomistic stress as both of them essentially measures the thermal lability of the peptide fragments.
Peptide fragments in the 318-323-328 K temperature segment (bottom panel, Figure 6) are intimately associated with the BT since a stress-temperature plot of each of these fragments (SI, Figures S1-S3), like peptides 248-255 ( Figure 5), can be associated with T BT . These peptide fragments can be collectively classified into two major categories. One that manifests reduced dynamic mobility after the BT. The peptides under this category constitute several long entangled loops centered around the capsid threefold axis (e.g. peptides 215-246). The other that manifests high dynamic mobility leading to BT. The peptides here are mainly b-strand that belongs b-sandwich fold of each capsid subunit from which the capsid spike protrude. A closer look further reveals that for the peptides 248-255 and 178-195 retains their dynamical mobility (cf. < 5 slope to fall from the maximum stress level) after breathing transition. In fact, peptides 248-255 are seen to hold onto the minimum slope in stress-temperature profile even up to 333 K, well above the discovered T BT (see Figure 5(A), bottom panel). Thus they can be identified as true hotspot region. They are the locations on the capsid structure from where structural disintegration initiate at T BT . As a result, water permeation rate through the capsid wall ( Figure 3) and so also is capsid volume (Figure 2) increases beyond T BT . The spatial range of the effect of these hotspot residues extend to 40 Å surrounding the fivefold axis, centering the VNP pore through which externalization of Nts begin. This result is again strikingly similar to HDX-MS experimental observation (van de Waterbeemd et al., 2017).

Simulation protocol
All MD simulations were performed in the GROMACS-4.5.5 suite (Hess et al., 2008) using the Amber ff99SB-ILDN force field with improved side chain torsion potentials (Lindorff-Larsen et al., 2010). Structure of the parvovirus MVM particle was taken from VIPER databank (Carrillo-Tripp et al., 2009) (pdb code 1Z14) (Kontou et al., 2005). The VP2 VNP consists of 60 identical copies (subunit) of protein chains. The capsid structure of the VNP is put into a cubic simulation box of dimension 32.3 nm Â 32.3 nm Â 32.3 nm and 828,768 TIP3P water molecules (Jorgensen et al., 1983) were added. Adequate number of Na þ ions are added to attain electroneutrality. Further, 100 mM NaCl is added, which roughly resembles the physiological ion concentration. No divalent cations are added in order to minimize the effect on stability due to them (Draper, 2004). The calculations were performed at neutral pH and the same HIS tautomer was used all over. Finally the simulating system is comprised of 2,998,794 atoms. The total thickness of the water reservoir around the VNP solute in all three directions was enough to ensure a bulk-like water behavior and smallest distance of the protein atoms from the bounding box is at least 10 Å. Long-range electrostatics were handeled with particle mesh Ewald method (Essmann et al., 1995), using a grid spacing of 1.2 Å and fourth-order interpolation with a 12 Å cutoff. For the short-range electrostatics and van der Waals interactions, a cutoff of 14 Å was applied. The integration time step was 2 fs with the LINCS method being used (Hess et al., 1997;Ryckaert et al., 1977). In order to check the accuracy of the chosen time step, for certain runs, 1 fs integration time step was chosen while we have used SHAKE algorithm. Note that for computing the mean-square displacement of H-atoms, the constraint algorithms were switched off. The nonbonded interactions are evaluated every 2 fs and electrostatics are updated every 4 fs.
The MD simulations were conducted in three stages: energy minimization, system equilibration and finally the production runs. The energy minimization of all the simulations were conducted with two algorithms; steepest-descent and conjugate gradient to reach a maximum force smaller than 300 kJ/(mol Â nm). The equilibration simulations were conducted in three stages; position restraint relaxation coupled to a heat bath, and NVT equilibration followed by NPT equilibration, each of 10 ns duration. The equilibrated structures, at temperature ranging from 100 to 440 K (from 100 K to 160 K at 10 K interval, thereafter up to 270 K at 5 K interval, 273 K, 298 K, 310 K, 318 K, 323 K, 328 K and 333 K, thereafter from 340 K to 440 K at 10 K temperature interval) were further used for subsequent 10 ns production runs under constant NPT conditions. The choice of the temperature range that started below the water freezing point and ended up above the boiling point ensured that the two transition temperatures can be captured with sufficient data available at both the temperature ends. The pressure of the systems was kept at an average of 1 bar using the isotropic Berendsen barostat (Berendsen et al., 1984) with a time constant of 2 ps and a compressibility of 4.5 Â 10 À5 bar. The system was coupled to an external heat bath with the velocityrescaling thermostat (Bussi et al., 2007). This thermostat however, suffers from 'flying ice cube effect', an artifact first discovered by Cheatham (Harvey et al., 1998). However, in GROAMCS suite, the problem associated with this energy partitioning artifact can be substantially reduced by the use of very long coupling times (we use a time constant of 0.2 ps) and frequent (after every 20 fs time steps) removal of the center of mass and global rotational motion of the system (using comm_mode ¼ angular keyword) (Harvey et al., 1998). To assess the effect of protein assembly on the onset of DT, we have also conducted NPT simulations under identical conditions for the monomer subunit (totalling 202,451 atoms) and the five protein subunits (totalling 928,740 atoms) surrounding a VNP pore. The convergence of each of simulations was checked by calculating the RMSD of C a atoms. In order to calculate water permeation through the semipermeable VNP (see below), simulations at temperatures, 273 K, 298 K, 310 K, 318 K, 323 K, 328 K and 333 K are continued for 200 ns. In order to assess the flexibility of the systems, the RMSD vs. time plots for the backbone atoms are prepared (SI, Figure S4). Note that due to the massively large system size (% 3 million atoms) and limited computational resources, we could not afford 200 ns long simulations at other temperatures than the seven temperatures reported here (cf. Figure 3). However, linear regression of the computed water permeation data and the exact tally of the discovered T BT value with that previously obtained from capsid's internal volume variations (cf. Figure 2), gave us confidence on the chosen few data points in Figure 3.
We have also conducted NVE MD simulations of the VNP to calculate the equilibrium atomistic mechanical stresses (see below) along the peptide fragments. In short, the varieties of simulations decipher the thermal stability, water permeability, and mechanical properties of the MVM capsid at atomic resolution.

Mean-square displacement (MSD)
Hydrogen being the lightest element will be the first to respond to temperature enhancement from cryogenic conditions and hence would be the first to exhibit the onset temperature of anharmonicity in atomistic displacement (Gupta et al., 2016;Pathak & Bandyopadhyay, 2019;Yoon et al., 2014). In view of this the MSD of protein H atoms (amide and all hydrogens) are calculated by, wherer i is the coordinate of the ith H atom, N is total number of the respective atoms, T is the length of trajectory used to calculate MSD, and dT is the MD time interval. This operational expression for MSD allows it to be estimated at various time windows of the entire trajectory, such that, its invariance with respect to time can be checked.

VNP volume
The capsid internal volume at varying temperature is determined by Richards' rolling probe method (Richards, 1977;Voss et al., 2006;Voss & Gerstein, 2010). The volume of the enclosed surface is determined by the discrete volume method of a rolling probe using 3 D grid of voxels. In brief, the solvent-accessible volume of the macromolecule is determined by reading the coordinates of its atoms and assigning to each atom a radius which is its VDW radius plus the radius of a suitable probe. The coordinates of the capsid forming atoms are placed on a cubic grid. The internal volume of the capsid is then determined by the set of all grid points that fall within the augmented sphere of any VNP atom. Finally, the volume is determined by multiplying the total number of occupied grid points by its voxel volume. The calculation time is linear in number of atoms, increases cubically with the probe radius and the inverse cube of the 3 D grid resolution. For all calculations we have used a probe radius 3 Å and gride size 0.4 Å (i.e. 16 voxels/Å 3 ). Note that the working expression for MSD (Equation (1)) ensures the calculated values are invariant to the trajectory lengths. However, VNP volume calculations are susceptible to the length of trajectory. In order to validate the discovered T BT value ( Figure 2) through VNP volume calculations, we have divided the 200 ns long trajectories at seven temperature points into 20 sub trajectories, each of 10 ns duration. Subsequently, VNP volume calculations are performed for 20 trajectory windows at each temperature and an error estimation is performed. The results are presented in SI Figure S5. As can be seen there fit to the data validates the discovered breathing transition at 318 K.

Water permeation
The VNP works as semipermeable membrane for water molecules and ions (Andoh et al., 2014;Perilla & Schulten, 2017). In order to calculate the water permeation rate one must define the boundary that divides inner and outer region of the capsid. Richards' rolling probe method is further utilized in this regard to find the boundary of the capsid. Namely, it can be found by the set of all unoccupied grid points that are next to the points defining the internal volume (see above). The surface boundary thus found reflects the atom arrangements of the VNP forming proteins along its wall. Exchange rates of water molecules in and out of the capsid is calculated by counting the number of water molecules crossing the boundary, b by, where i runs over the waters in a given side of the capsid (interior or exterior) during t À Dt to t: We choose Dt ¼ 1 ns, the time between trajectory frames of 200 ns long simulation at each temperature. The reference starting time-frame was also varied to check the convergence of permeation events. Slope of number of water molecules going in and out of the boundary versus time provide the rate of permeation.

Atomistic mechanical stress
The calculation of local stress field in an atomistic simulation of capsid is essential in order to both characterize structural disorder as the temperature of the solvent bath increases and to interpret the results in light of the 'breathing motions'. Here we adopt an original formulation for the atomistic stress calculation based upon the idea of a volumetric partitioning of the bulk stress tensor (Basinski et al., 1971), which have been elegantly applied in several other occasions (Fenley et al., 2014;Hatch & Debenedetti, 2012;Zimmerman et al., 2004). In this formulation, the basic assumption is that the continuum theory of elasticity continues to characterize the stress in a small volume element X a about atom a even when stress distribution is inhomogeneous as is expected in the polyhedral atomic assemblage of a capsid. Thus the microscopic stresses can be found at various length scales varying from atoms to individual residues, to peptide fragments and to the whole capsid. In particular, if V is interatomic potential in an additive force field setup, the forcef ab exerted by atom b upon atom a is given as f ab ¼ À oV=or ab À Ár ab =r ab then atomic stress tensor can be written as, where m a is the mass andṽ a is the velocity of atom a and the interatomic separation,r ab ¼r a -r b : The first term in Equation (3) accounts for momentum transfer by translation of atoms through space, which in the local atomic frame of coordinate system would disappear. Hence in the local (Lagrangian) frame of atoms, Equation (3) can be simplified to contain only the second interatomic force term. The multicomponent tensor containing the interatomic force term at each atom in the capsid is evaluated by averaging the principal stresses, which is simply one-third of the trace of stress tensor (Hatch & Debenedetti, 2012). In order to calculate the atomistic mechanical stress of the capsid in explicit water solvent, NVE simulations of at least 50 ns (1 ps time step) are performed using starting structures extracted from the previously conducted NPT simulations at various temperatures. The convergence of the NVE equilibrium simulations are checked by computing the root mean square deviations of the C a atoms of the VNP protein residues from their initial starting structures. Each equilibrated structures are further subjected to 5 ps NVE simulations with 0.1 fs time steps. Finally, the total stress accumulated over 1 ps time window in each of the capsid residues are calculated. The comparison of accumulated stresses in the five time windows are made in order to verify the mitigation of the stress fluctuations induced by thermal fluctuations. The time-averaged total mechanical stress of the full VNP structure (with M number of atoms) is found by r capsid ¼ 1 ðTÀdtÞ Ð TÀdt 0 P M a¼1 r a t ð Þdt: In summary, 10 ns NPT MD runs are conducted at temperatures ranging from 100 to 440 K (cf. Figures 1 and 2). Out of these, runs at 273 K, 298 K, 310 K, 318 K, 323 K, 328 K and 333 K are extended up to 200 ns (cf. Figure 3). NVE simulations (first for 50 ns with 1 ps time step and then for 5 ps with 0.1 fs time step) were also performed for the structures at 273 K, 298 K, 310 K, 318 K, 323 K, 328 K, 333 K, and 400 K temperature points (cf. Figures 4-6).

Conclusion
We have carried out large-scale MD simulations of a VNP to investigate the heat-induced transitions, namely, DT and BT. Although both the transitions result from the collective atomistic motion of the shell, the associated temperatures are widely separated from one another (T BT >T DT ) by $ 100 K. We found that at T DT , DT brings forth the flexibility in the shell, which ceased to exist at cryogenic conditions. Whereas at sufficiently higher temperature, T BT , BT brings it in a functional state, wherein global structural rearrangement of VNP proteins reaches a critical state such that release of N-termini through the pore can initiate. While for the shell, heat induced changes in the fluctuation of hydrogen atoms quantifies the onset temperature of DT, a host of macroscopic parameters, such as changes in internal volume, rate of water permeation through the semi permeable surface, and the elastic properties in response to temperature variation uniquely determine its T BT . The mechanical response of capsids (Carrasco et al., 2008;Castellanos et al., 2012;Kol et al., 2006;Medrano et al., 2019;Michel et al., 2006) to environmental insults, their stability, and the threshold of their thermal lability are major issues pertaining to bio-inspired engineering of nanocontainers (designed for drug delivery, for example). Although the macroscopic parameters monitored here are evaluated from all-atom MD simulations, they lack the atomistic details of the shell in the sense that they remained silent about the breathing motions of the shell in which the different capsid proteins fluctuate between different conformations and only a couple of them (hotspot regions, see SI Figure S6) orchestrate the global structural rearrangement necessary for the translocation of N-termini through the pore. To gain insights into the heat-induced breathing mechanism, atomistic mechanical stress-strain calculations are performed in the NVE ensemble by evaluating the corresponding virial of each atoms of the VNP. The results form atomistic stress calculation can be summarized as: (i) the distribution of stress on the VNP structure is highly inhomogeneous, and (ii) total (sum of tensile and compressive) stress-temperature profile of the complete VNP also exhibit BT at T BT intimately close to the value extracted from measured macroscopic properties. The last result re-establishes the close connection between macroscopic observables and microscopically atomistic stress. Further, based on the thermal response of the stress levels of various peptide fragments, they are classified into several categories that helps identify the hotspot regions. The discovered T BT value and the breathing mechanism are in remarkable unison with deuterium uptake (HDX-MS) experimental results since both of them in essence measures the thermal lability of the peptide fragments. In particular, the most thermo-labile peptide fragments (178-195 and 248-255) of the MVM particle are not found to be the pore forming residues, rather they are at least 40 Å spatially separated from the capsid pore. It would be interesting to investigate whether the most thermo-labile peptide fragments identified in this work will remain invariant with the variation of other environmental factors like, pH. Coarse Grain simulation methods could show a good approximation to the event. In view of the recent upsurge in interest in the synthesis of bio-nanocontainers for their potential applications in material/medical science, the present in silico study opens up an additional avenue to infer the structure-stability relationship of the assembled systems, at a level of detail, which is not always accessible to experimental methods.
The present work is carried out in nonpolarizable water media despite the fact that polarizability of water is different in the bulk and in protein. As a result, the observations made here on the studied phenomena will change in a more realistic polarizable medium, although their qualitative trend (i.e. T BT >T DT ) is expected to remain intact (Pathak & Bandyopadhyay, 2019). We hope to address this issue in the near future.