Temperature-mediated switching of protectant-denaturant behavior of trimethylamine-N-oxide and consequences on protein stability from a replica exchange molecular dynamics simulation study

Abstract The detailed mechanism of protein folding–unfolding processes with the aid of osmolytes has been a leading topic of discussion over many decades. We have used replica-exchange molecular dynamics simulation to propose the molecular mechanism of interaction of a 20-residue mini-protein with urea and trimethylamine N-oxide (TMAO) that act as denaturing and protecting osmolyte, respectively, in binary osmolyte solutions. Urea is found to exert its action by interacting directly with the protein residues. Temperature tolerance of TMAO’s action is particularly emphasised in this study. At lower range of temperature, TMAO acts as a successful protein protectant. Interestingly, the study discloses the tendency of TMAO molecules to prefer self-association at the protein surface at elevated temperature. A greater number of TMAO molecules in the protein hydration shell at higher temperature is also observed. Dihedral angle principal component analysis and free energy landscape plots sampled all possible conformations adopted by the protein that reveal highly folded behaviour of the protein in pure water and binary TMAO solutions and highly unfolded behaviour in presence of urea.


Introduction
The folding and unfolding of proteins is of vital importance for the regulation of their biological activities. The folding of proteins into their native states is the most fundamental example of biological self-assembly. Hence, understanding of this process will help to elucidate the process of evolutionary selection of properties of a biological system for functional advantage. Native state of a protein refers to the state of the highest thermodynamic stability under physiological conditions. If the folding process of a protein from its unfolded state undergoes all possible conformations, it will take an astronomical length of time to reach to the native state. However, literature [1,2] revealed that folding process does not involve a series of mandatory steps between specific steps, but rather a stochastic search path is followed. The compactness of the three-dimensional structure of a protein depends intrinsically on the sequence of amino acid residues as well as on the multiple contributions from the cellular environment. The uniqueness of energy landscape of a particular protein also pertains to amino acid sequence and this leads natural selection to evolve so that the protein gets native folded state rapidly and efficiently.
The challenges in the conformational sampling of complex biological systems prohibit the complete exploration of energy surfaces. In literature, many significant efforts have been reported till date for developing efficient simulation methods to overcome the trapping (of the system of interest) at low-energy minima which may prevent sampling of some important regions of the energy landscape. Though  Exchange molecular Dynamics (REMD) in 1999 by Sugita and Okamoto [3], REMD has been a very popular approach to carry out biomolecular simulations, its application to large systems still remains to be exploited.
Protein is sensitive to the surrounding environmental conditions such as temperature, pressure and assorted molecules comprising the volume around it [4]. Some small molecules that accumulate in high concentration inside the cell in response to osmotic stress are known as osmolytes. Nature exemplifies certain evidences of using protecting and denaturing osmolytes to restore cellular functions. For example, certain marine animals have adapted life at very high pressure and salinity using osmolytes having denaturing effect such as urea [5] and protecting osmolyte such as betaine, trimethylamine N-oxide (TMAO), etc., to counteract the effect of denaturant. [6] Urea is an extensively studied denaturing osmolyte, which has been postulated to exert its effect either directly [7] (by binding to the protein) or indirectly [8] (by altering the solvent structure which in turn offers denaturing effect on the protein). Some studies attributed the urea denaturation to both "direct" and "indirect" mechanisms [9,10] while others emphasised the role of hydrogen bonding between urea and protein. [11,12] However, the exact mechanism of urea-induced denaturation is not adequately established till date. [13] The osmolyte TMAO is found in the urea-rich cells of elasmobranchs and coelacanths to protect the protein against deleterious effect of urea. TMAO is considered as counteracting osmolyte which imparts its effect on both structural stability and functions of the protein.
Studies based on TMAO-water interaction showed that the unfavourable interaction of TMAO with protein backbone is dominantly responsible for its stabilisation effect. [14] Further, studies on the molecular interactions between TMAO and protein functional groups also give evidence of some favourable and unfavourable contributions towards the overall stabilising effect of the osmolyte. [14,15] TMAO also may act as a crowding agent leading to entropic stabilisation mechanism of protein through excluded volume effects. [16] Here we note that a recent study looked at the origin of the contrasting behaviours of the aqueous solutions of isosteric molecules TMAO and TBA (tertiary butyl alcohol) on the protein stability. [17] It has been argued that in aqueous TMAO solution the stability of the native conformation of the protein arises because of the fact that the addition of TMAO to pure water increases density of the solution which in turn increases magnitude of the solvent-excluded volume effect. On the other hand, aqueous TBA solution has density lower than that of pure water and this leads to a reduction in the magnitude of the solvent-excluded volume effect. A much larger volume packing density of ternary aqueous solution of TMAO and urea than that of pure water, which makes the cavity creation process much more costlier for the former, has also been considered as a prime factor for the stability of the native conformation of the protein in ternary mixed TMAO and urea solution. [18] Simulation studies on TMAO-water mixture provide strong support of the presence of compactly bound water molecules around the polar head group of TMAO molecule. [19] Though a lot of research works have been devoted, still study on this crucial osmolyte is captivating new area of surveillance everyday. In literature many experimental evidences also exist showing Urea-induced denaturation, TMAO-induced protection and also counteraction of TMAO against urea action. [20,21] These findings are substrate specific based on intermolecular interaction mechanism. These studies inspired us to investigate and characterise the protection and denaturation mechanisms of osmolytes TMAO and urea, respectively, at molecular level using REMD simulation technique. To bring out the detailed mechanism of protein-solvent interactions and their roles in overall conformational behaviour of the protein are the core objectives of this study.
The rest of the paper is divided into three sections. Section 2 gives a description of models and simulation details, Section 3 presents a discussion of the results and conclusions are summarised in Section 4.

Models and simulation method
One 20 residue mini-protein Trp cage ( Figure 1) having amino acid sequence NLYIQWLKDGGPSSGRPPPS (PDB code 1L2Y) was chosen for the study. This mini protein has two positively charged side chains (Lys 8 and Arg 16) and one negatively charged side chain (Asp 9). Hence one chloride ion was added to neutralise the system. The protein was solvated with TIP3P water [22] in a cubic box and with respective number of cosolute molecules for preparing the different systems. To study the denaturing effect 5.8M urea was used, whereas 2.9M TMAO was used as protecting osmolyte. The details of simulation of all the systems simulated are presented in Table 1. Respective parameters for urea model, [23] Kast model of TMAO [24] and TIP3P water are included in Table 2. For Trp Cage, AMBERff94 force field parameters were used. Here we would like to mention that a separate simulation run was performed in which the Kast model of TMAO was replaced by the equal number of Osmotic model of TMAO molecules. [13] The consideration of two different models for TMAO molecules would allow us to verify the model dependency of our results if any.
For each of the simulation run the system was subjected to 2500 steps of SD (steepest descent) minimisation and was subsequently equilibrated at constant pressure (1 atmosphere) and temperature (300K) for 300 ps. This results into a system of adequate density with corresponding edge lengths. Further, an equilibration at constant volume was performed for 200 ps. In order to set up the distribution of temperatures for replicas, the system resulted from the above mentioned steps were equilibrated further at seven different temperatures (viz. 275, 325, 375, 425, 475, 525 and 575K). The average energy values obtained from these simulations were then fitted with a polynomial equation which was solved repeatedly to find out temperature distributions such that the acceptance probability was 20%.
AMBER12 [25] molecular dynamics (MD) package was used to perform simulation of the replicas which were first equilibrated separately for 50 ps at constant volume at their respective temperatures. All the systems were then subjected to 80 ns simulation run each. Exchanges were attempted every 500 integration steps. After the initial equilibration period, coordinates were saved for every 2 ps over the whole production simulation time (in NVT ensembles). All bonds involving hydrogen were constrained with SHAKE [26] with a 0.00001 Å tolerance. Temperature was controlled by Langevin dynamics [27] with a collision frequency of 1 ps −1 . A 9.0 Å cut-off was used for all non-bonded interactions. The electrostatic interactions were treated by particle mesh Ewald summation method. [28] Periodic boundary conditions were implemented in all three directions which render the simulation of fully solvated central cell in the environment produced by the repetition of the cell in all directions. The replica-exchange molecular dynamics (REMD) algorithm of Sugita and Okamoto [3] arises by the application of the parallel tempering method to MD simulation, where multiple copies (or replicas) of identical systems are simulated simultaneously and independently at different temperatures for a chosen period of time or MD steps. Then, periodically stateexchange moves are attempted, where two neighbouring replicas exchange their thermodynamic state (i.e. temperature). The exchange processes among the replicas are designed to maintain Boltzmann distribution. The acceptance rule for each stateexchange moves between two neighbouring states i and j is chosen to be where E i and E j represent the configurational energy of the system of replica i and j, respectively. The temperature spacing between the replicas is chosen such that sufficient overlapping of the energy distribution between neighbouring replicas is allowed. If the exchange is accepted, the replicas will swap their temperatures and velocities will be scaled accordingly. On the other hand, if the exchange is rejected, no swapping of temperature takes place and each replica continues with its existing trajectory.
Trajectories obtained from REMD simulations were analysed using both AMBER ptraj module as well as VMD (Visual MD package) [29] to study the effect of cosolvents in the folding/ unfolding equilibrium of Trp cage. For the sake of promoting an equilibrium view, analyses were done on the trajectory obtained after achievement of the proper convergence.

Checking of convergence
In order to check simulation convergence, helicity of all the residues are calculated for every time step of a particular temperature, viz. 284K of all the three systems. We have used VMD [29] for the calculation of backbone dihedral angles phi( ) and psi( ) and for calculating helicity we set the criteria as −90 • ≤ ≤ −36 • and −72 • ≤ ≤ −0 • . The plot of helicity for first 60 ns is shown in Figure 2. Homogeneity of the helical behaviour throughout the timescale gives evidence of the convergence for respective system. It is observed that after first 20 ns, each of the system reaches equilibrium.

Probability distribution of radius of gyration
To quantify the degree of folding of the protein, radius of gyration (R g ) is computed. R g measures the effective size of the protein chain. In order to compare the folding behaviour of the protein, probability distributions of radius of gyration of C α atoms of Trp cage for all the systems at low-(284 K) and high-temperature (≥590 K, see Table 1) limits are plotted (see Figure 3). Note that, since our results for 284 K temperature is very similar to that of at 300 K temperature (not shown) and the former temperature (i.e. 284 K) is common to all the systems, for low-temperature limit we show the results for 284 K. It is observed that pure water and TMAO-containing systems lead to a narrower distribution of R g values at low temperature (when compared to that of urea-containing system) indicating the decrease in configurational space. This suggests contributing effect of TMAO at low temperature which helps the protein to remain in compact state, i.e. folded conformation. The R g value of the protein for the system containing Osmotic model of TMAO shows a very similar behaviour as that of the system containing Kast model of TMAO (see Figure 1 in Supporting Information). But for the system containing urea, the protein shows both folding and unfolding behaviours (since R g values spread over a wide range) due to denaturing contribution of urea at both the temperature limits. At high temperature, ( Figure  3(b)) a much wider distributions of R g values indicate unfolded conformations for all the systems.

Melting curves
In order to obtain the temperature-dependent behaviour of the protein, it is rational to examine the dependence of native fractions on temperature. According to Anfinsen [30], a particular conformation of a protein assumes the state which is thermodynamically most stable under any set of specific conditions. In our study, we use a threshold of RMSD ≤ 2.3 Å to figure out folded state. This criterion is set from the distribution of RMSD values of the α-helical part of the protein at the lowest temperature (not shown here). RMSD is calculated according to the following equation: is the position of ith atom at time t and m i is the mass of the ith atom.
The folded fraction of the protein in presence of water and different milieu as a function of temperature is shown in Figure 4. It can be seen that with increasing temperature, folded fraction decreases for each of the systems. For pure water system, folded fraction does not change till 400 K and melting temperature approaches 500 K. Here we note that the experimental melting temperature of the protein Trp cage is about 315 K which is confirmed by various experimental methods like NMR, CD, DSC, fluorescence and UV-Raman by a number of research groups. [31][32][33] However, in simulation studies the melting temperature of the protein Trp cage varies from 360 K [34] to around 450 K [35,36] depending upon the force field parameters used.
Different shapes of the curves that exhibit the melting behaviour of Trp cage are affected differently depending on the nature of the cosolvent present. Addition of urea decreases the melting temperature to below 400 K. In presence of TMAO, folded fraction decreases gradually with increasing temperature and melting temperature goes to around 450 K. In Supporting Information Figure 2(a), we have compared melting curves of binary TMAO systems with both Kast and Osmotic models of TMAO. The results confirm that the action of the two different models of TMAO on the protein are very similar.
Over the entire range of temperatures, the most stable behaviour of the protein is observed when there is no cosolvent. Since TMAO is known to act as protein stabiliser, hence the most stable behaviour of the protein is expected with 2.9M TMAO solutions. But, discussions of radial distribution functions (rdfs) and hydrogen bonds (discussed later) give evidences of highly interactive nature of TMAO with protein over the whole temperature range. In the binary system containing 5.8M urea, the behaviour of the protein is as per the expectation. From very low temperature itself unfolding starts and it continues up to high temperature. Figure 4(b) shows the respective values of change in Gibbs free energy, G u , of unfolding for all the systems studied against temperature. G u can be calculated as: where K is the equilibrium constant and x f is the folded fraction of the protein. Insertion of TMAO in the water-protein system causes an increment in the stabilisation of the protein than in presence of urea with a corresponding increase in free energy values. For example, the free energy of unfolding G u = 5 kJ/mol at 284 K for binary TMAO system which approaches to 0 kJ/mol for binary urea system. At this temperature, the maximum stable form of the protein is retained in absence of any cosolvent with free energy of unfolding value of 13 kJ/mol. Here we note that a comparison between the G u values for two different models of TMAO reveals that at lower limits of temperature, in the system containing Kast model of TMAO the folded fraction is slightly higher than that for Osmotic model of TMAO (see Figure 2 (b) of Supporting Information section).

Free energy landscape
The free energy landscapes (FELs) are plotted to obtain a detailed comparison of energy surfaces of all the systems. For this we have used the following equation: where P(V) is the probability distribution of the respective coordinate (V) obtained from the simulation trajectory. P max denotes the maximum of the distribution, which is subtracted from P(V) to ensure that G = 0 for the lowest free energy minimum. In our study, G is plotted vs. C α RMSD for respective temperatures. The different shapes and clearly distinguishable states in FELs are indicative of different intermediates throughout the process of folding to unfolding. FEL plotted ( Figure 5(a)) for system comprising of protein and pure water shows a distinctive region of folded state which is stable up to ∼ 480 K. Addition of urea leads the FEL ( Figure  5(b)) to take a diversified shape containing both folded basin and unfolded basin. In presence of TMAO, FEL shows narrower distribution of RMSDs throughout the whole temperature range ( Figure 5(c)) than the urea-containing system. A comparison of all FELs clearly indicates that the stabilisation of the protein (with RMSD ≤ 2.3 Å) is maximum in absence of any cosolvent. FEL for the binary TMAO system with Osmotic model is also shown in Supporting Information Figure 3  conformational space for system consisting of Kast model of TMAO is found to be narrower than that of Osmotic model.

Dihedral angle principal component analysis
For principal component analysis (PCA) of a system containing M atoms, the correlated internal motions are represented by the covariance matrix: To obtain Dihedral angle principal component analysis (dPCA) [37] we have performed PCA on sin-and costransformed dihedral angles. The free energy profile is plotted according to Equation (4) which is highly rugged with multiple number of minima.
Three different free energy surfaces (see Figure 6) at 284 K temperature are plotted using 'Carma: a MD analysis program' [38]. The figure shows free energy as a function of first two dihedral angle principal components. dPCA can clearly resolve numerous minima of the free energy surface. The structures of the protein obtained from the clusters of the corresponding minima are also presented. The different landscapes efficiently resolve the conformational heterogeneity of Trp cage under different environments.
In Figure 6(a), the five minima belong to folded conformations in which the protein spends ∼ 66.3% of the entire simulation time. Similarly 8 and 4 number of minima are found in Figure 6(b) and (c) with ∼ 50.8% and ∼ 59% occupation times, respectively. High-and low-energy barriers (with different energy values) separate the basins. In Figure 6(b), the minima ranges from folded to partially unfolded one. The minimum corresponding to unfolded conformation shows high energetic barriers. Figure 6(c) shows four distinct minima with folded conformations which are well separated from each other. Lowenergy barriers in this figure depict easy transitions between the closely related conformations.

Solvation of protein
Water plays a significant role in determining the structural stability and the function of protein. [39,40] Presence of any large macromolecular solute affects immensely on both structural rearrangement and mobility of the solvent. Literatures in the field of protein crystallography and MD simulation aid to confirm this effect. [41][42][43][44][45][46] Global picture of protein hydration [47] depicts four distinct regions around it: (a) the solvent-free region inside the protein where there is no solvent molecules; (b) surrounding the solvent-free region there is a transitional second region of inter-penetration of protein and solvent; (c) this is followed by a region of local maxima of solvent density and (d) a fourth region is the remote region with no distinct density maximum. Computed site-site radial distribution functions (rdfs that depict the solvation pattern of a reference molecule) of oxygen atom of water around the protein heavy atoms for all the systems are shown in Figure 7. The closest peak appears at about 2.8 Å with a small peak height, which is an indication of interaction of water and protein. The second peak appears at about 3.8 Å with considerably higher peak height than the first one portraying non-hydrogen bonding interaction. A third peak with small peak height emerges at about 4.8 Å demonstrating isolated hydration sites around the protein (with a much weaker interaction).
At high temperature (≥590, see Table 1), the first peak disappears and the second peak height increases for all the systems (see Figure 7(b)). This suggests that the hydrogen bonding interactions (between protein and water) become less significant though water density and local hydration sites in the first solvation shell (FSS) of protein increases. But both at high and low temperatures, water density in close proximity of the protein is lower than that in bulk which narrates the exclusion of water molecules from the protein surface. This is in consistent with the fact of natural tendency of protein to remain in the folded state as long as water solvation persists. [48] Difference in hydration behaviour of the protein in presence of different osmolytes is also quite evident. In binary urea-water system, at low temperature the first peak height is raised, i.e. more number of water molecules are incorporated around the protein through hydrogen bond interaction and diminished peak heights of second and third peaks reveal that non-hydrogen bonded water number decreases indicating that those proximities are filled by urea molecules. In binary TMAO system, all the peak heights increase to a little extent portraying enhanced water solvation due to the presence of TMAO.
Calculation of the number of water molecules around the different residues of protein and the average number of proteinwater hydrogen bonds (see Tables 3 and 4) also reveal some facts about protein hydration behaviour. Water molecule having its oxygen atom within a distance of 3.5 Å from any heavy atom of the protein is considered to be inside the FSS. And a maximum D-A distance of 3.5 Å and a maximum D-H-A angle of 45 • (following earlier work [48]) are used as criteria to find the number of different types of hydrogen bonds (where D and A represent donor and acceptor atoms, respectively).
It is observed that there are more number of water molecules in the FSS of the protein at high temperature than that of low temperature for each of the systems, implying more exposure of the protein to solvent water at high temperature (see Table 3). Calculated numbers of heavy atoms of urea and TMAO in the FSS of protein for respective systems are also included in Table 3. When osmolyte is present in the system, the number of water molecules inside the FSS decreases. Further, these data give affirmation regarding more tendency of urea molecules to remain in the close proximity of the protein than TMAO molecules. With increasing temperature, urea molecules disperse away from the FSS while slightly more number of TMAO molecules tend to accumulate inside the FSS. It is also noticed that urea molecules tend to interact with the protein mostly through its nitrogen and oxygen atoms, while TMAO prefers to interact through its methyl carbon atoms. These observations are in accordance with the previous results reported elsewhere. [49]  For pure water system, the increase in number of water molecules in the FSS of the protein with increasing temperature clearly indicates the denaturation of the protein, which leads to extended structure and hence more number of water molecules are accessed inside the FSS. When TMAO is added to the system, a very few molecules of TMAO also come inside FSS along with water molecules. At high-temperature limit, the number of both water and TMAO molecules increases indicating extended structure of the protein. For the system-containing urea, it is found that less number of water molecules are present in

System
Temperature the FSS (when compared to that for pure water system) due to the penetration of urea molecules inside the FSS at both temperature limits. Interestingly, the increase in temperature shows two opposite effects in regard to the number of each of the osmolytes in FSS. In specific, the number of urea molecules in the FSS decreases and that of TMAO increases as temperature is increased. The calculation of percentage of water molecules present in the FSS of the protein, which are capable of forming hydrogen bonds with the protein, shows that (see parentheses of Table 3) in presence of urea the water molecules are slightly more susceptible to form hydrogen bonds with the protein when compared to pure water system. And in presence or absence of TMAO the behaviours of water molecules to form hydrogen bond with the protein are similar. As expected, for all the systems, hydrogen bonding capacity of water within FSS decreases with increasing temperature.
Analyses of average number of hydrogen bonds (as shown in Table 4) formed between Trp cage and water show that in absence of any cosolvent the protein is least likely to form hydrogen bond with the water molecules (when compared to osmolyte-containing systems). This proves that water does not act as a good solvent to solubilise protein. It helps to elucidate the idea of spontaneous folding of protein in aqueous environment. [48] Since urea molecules can form considerable number of hydrogen bonds with protein they deliver denaturing effect onto the protein. Average number of hydrogen bonds formed by the protein with water is also found to be increased in presence of urea. As a result there is an overall increase in the hydrogen bond numbers that the protein engages in aqueous urea solution when compared to that for pure water system. [49] At high temperature, the number of hydrogen bonds formed by urea with the protein is decreased to half while the number of protein-water hydrogen bonds remains close to that of the low-temperature value. The protein rarely forms any hydrogen bond with TMAO. But the number of protein-water hydro-  gen bond increases in presence of TMAO. Nevertheless, for binary systems the total number of hydrogen bonds formed by the protein with water and osmolyte (i.e. either urea or TMAO) is much higher for urea containing system. To obtain a peer view of TMAO's action on the protein, hydrophilic and hydrophobic residues are considered separately to study their interactions with the active sites of TMAO (see Figure 8). The central atom of TMAO (i.e. nitrogen) keeps the maximum distance from the protein residues. For binary TMAO system, carbon atoms of TMAO prefer hydrophobic sites of the protein and its oxygen atom prefers hydrophilic sites. These distribution functions further reveal the accumulation of more number of TMAO molecules near to the hydrophobic sites of the protein at higher temperature whereas increased temperature has little effect in TMAO accumulation around the hydrophilic sites.
These imply that TMAO's preferential accumulation near to the hydrophobic sites of the protein.
In this context spatial distribution functions involving TMAO methyl carbon and the protein are plotted (see Figure 9). The figure describes low occupancy of TMAO methyl carbon atomic sites near the surface of the protein at low temperature. However, at high temperature, the tendency of TMAO molecules is raised to occupy a nearby position to the protein.

Water structure arrangement
Since TMAO does not interact directly with the protein particularly at low temperature, it becomes important to investigate the behaviour of water in presence of TMAO. Literatures show that TMAO creates an environment where unfolded state of the protein is undesirable in comparison to the folded state due to solvophobic effect of TMAO on the protein. [50,51] In view of this, it is important to examine the role of added osmolytes on the water structure through distribution functions involving water molecules. Water oxygen-oxygen rdf ( Figure  10) at low-temperature limit for each of the systems shows a sharp peak at 2.8 Å which is a characteristic of highly populated hydrogen bonded first neighbour. In comparison to this, the second and third peak height diminishes indicating that tetrahedrally located second neighbour and third neighbour are less pronounced. Presence of TMAO is found to be contributive for water molecules to become exceedingly hydrogen bonded along with well-heeled second and third neighbours. This makes the water network more structured which in turn makes the protein backbone less susceptible for hydration, i.e. leaving the backbone intact with its folded form. Rdf at high-temperature limit shows that the first peak appears at the same position due to hydrogen bonded first neighbour but second peak appears at a different position (about 5.6 Å). Further, the height of the first peak decreases as temperature is increased suggesting less number of water-water hydrogen bonds.

Self association of osmolytes
With increasing temperature if the molecules tend to stay apart from each other it will result into decreased peak height of the first peak of the rdf of central atom in comparison to the peak height at low temperature. But TMAO is found to exhibit exceptional behaviour regarding this aspect. Study of distribution function of TMAO-nitrogen around TMAO-nitrogen (as displayed in Figure 11(a)) shows that the first peak appears with magnified peak height at 5.1 Å for high-temperature ensembles indicating that TMAO molecules prefer to be in close proximity of other TMAO molecules at high temperature than that at low temperature. However, for the urea-containing system, Figure  11(b) shows a higher tendency of urea molecules to remain more associated with other urea molecules at low temperature with an elevated peak at 4.4 Å. At high temperature, this peak height diminishes revealing decreased tendency of self association of urea molecules. This illustrates the behaviour of urea molecule to exert its denaturing effect uninterruptedly at high temperature. Further, in order to confirm the greater existence of TMAO cluster at higher temperature we have estimated mean cluster size at the surface of the protein at the two temperature limits. The last 30 ns of our simulation run is used for this study. After achievement of convergence for the binary TMAO solutions, clustering behaviour of TMAO molecules is studied (at interval of 5 ns). TMAO molecules having any of its heavy atoms within 3.5 Å of the protein are selected and among those molecules  the one for which methyl-methyl (of TMAO) site distance falls within a cut-off distance of ≤4.5 Å are considered to form a dimer. Here, we note that we have not observed any TMAO cluster of sizes greater than two. From their number of occurrences, the probability of formation of binary cluster is calculated. The mean cluster size can be calculated as [52]: where n TMAO is the number of TMAO molecules in a given cluster and P TMAO is its probability of formation. Figure 11(c) is the plot of n TMAO as a function of time scale at the two temperature limits. The figure clearly depicts high occurrence of TMAO binary clusters at higher temperature than at low temperature at the protein surface. In this context, it is worth mentioning that room temperature MD simulation study of aqueous TMAO solution showed the non-existence of TMAO self-aggregation. [19]

Preferential local distribution
The different degree of solvation of the protein in presence of added osmolytes is characterised by estimating the relative local distribution of water and osmolytes around the protein. Following previous works, [53,54] the time-averaged normalised ratio (g i ) which provides information about how protein interacts preferentially with one type of cosolvent relative to another is calculated by: where n i and N i represent the total number of species 'i' defined in the local domain of radius 'r' of the protein and the total number of the species 'i' in the entire simulation box, respectively. When the vicinity of the protein is richer with component 'i' the value of the normalised ratio g i goes above 1 and the value below 1 corresponds to the depletion of the component in the nearby region of the protein. A comparative survey at low-and high-temperature limits (as shown in Figure  12) shows that urea prefers to stay nearby the protein and TMAO shows preferential exclusion from the protein surface at low temperature. For binary urea solutions, g U C exceeds value 1 up to a large distance from protein surface and then it starts to decrease. This indicates that urea molecules prefers to interact with the protein surface. At high temperature g U C is found to exceed 1 slightly since solvation shell of the protein becomes heavily loaded with water (as is evident from the discussions of Table 3). At low temperature g T N does not exceed 1, however at high temperature it increases slightly. Previous discussions of solvent count (Table 3) also fit with these findings. Therefore, changing behaviour of TMAO (with temperature) is found to play a dictating role in determining folding equilibrium of the protein.
In Table 5, we compare the relative population of water in the vicinity of the protein in comparison to the bulk. To fulfil this approach the ratio of water:osmolyte is calculated for all the systems near the surface of the protein (within a distance of 3.5 Å i.e. FSS) as well as that of the bulk for both the temperature limits. The maximum relative population of water on the protein surface is obtained in presence of TMAO at low temperature. This observation is in consistent with the fact of preferential exclusion of TMAO from the protein surface. [48] At higher temperature the decreased ratio of water:TMAO is due to accumulation of TMAO molecules at the protein surface. For the urea-containing system at low temperature, the relative population of water is low indicating preferential solvating nature of urea which prevents direct occupation of water molecules. But at high temperature limit, the ratio increases (unlike the case of TMAO) which is an indication of higher relative population of water molecules near the protein surface for urea-containing system. For both the cases, the relative populations of water molecules are found to be richer in the bulk region than the surface of the protein.

Summary and conclusions
We have carried out REMD simulations of Trp cage both in pure water and in presence of denaturing and protecting osmolytes. Different analyses clearly depict a stable and folded conformation of the protein at low temperature in pure water and in presence of TMAO. At high temperature, the protein becomes unfolded. But the degree of unfolding is different depending on the surrounding osmolyte, which is quite evident from the melting curves. Pure water system shows denaturation to small extent in comparison to the binary osmolyte systems. Urea causes maximum unfolding. Binary systems of urea and TMAO exhibit similar distributions of radius of gyration (shown in Figure 3(b)). In this context, the discussion from the calculations of rdf of central heavy atoms of TMAO ( Figure 11) is found to be helpful, which explains self-associative tendency developed by TMAO molecules at high temperature.
Calculations of solvent count and rdfs (at low temperature) reveal that water solvation is enhanced in presence of TMAO in the protein surface. The highly ordered structure of water network and highly stable state of protein in presence of TMAO are not observed in our study which are well established in other regular MD simulation studies. [50,51,55] Urea is found to exert its effect by direct interaction with the protein. Urea molecules possess high tendency to occupy solvation shell of the protein at low temperature. It is also found that in presence of urea, water molecules in FSS are very susceptible to form hydrogen bonds with the protein. But with the rise of temperature, number of urea molecules depletes from the solvation shell of the protein.
The difference in the action of urea and TMAO are discussed extensively by Berne and co-workers. [56,57] They suggested that both urea and TMAO interacts strongly with the polymer chain. The strength of van der Waals (vdW) interactions of a polymer chain with its surrounding osmolyte determines the effect of osmolyte on the polymer chain. Sufficient vdW interactions between polymer chain and TMAO disfavour extended conformation of the chain, while sufficiently strong interactions between polymer chain and urea favour extended conformation. It has also been suggested that thermodynamic binding affinity of TMAO molecules is higher for the collapsed conformation of a hydrophobic polymer than the extended conformation and the reverse is applicable for urea molecules. Effect of osmolytes on polymer conformation is determined by their accumulation or depletion in the FSS of the polymer. FELs of Trp cage for different systems signify that folding/unfolding equilibrium remains highly steady (up to 480K) in pure water. In binary TMAO solution, overall conformational space decreases than that of binary urea solutions. Protein protecting behaviour of TMAO at lower temperature is also reflected from FEL and dPCA plots. Highly rugged energy surfaces and representative structures from dPCA analysis for binary urea solution reflect unstable be-haviour of the protein in presence of urea. Furthermore, already existed osmotic model [13] of TMAO, where TMAO-TMAO Lennard-Jones interactions between different atomic sites are scaled down, is compared to the TMAO model employed in this study. We have found that Kast model is slightly more effective when compared to Osmotic model to provide protection to the protein when used in the same urea:TMAO ratio. These observations are in accordance with our previously reported results. [58]