Molecular Dynamics Simulations of the Temperature Induced Unfolding of Crambin Follow the Arrhenius Equation.

Molecular dynamics simulations have been used extensively to model the folding and unfolding of proteins. The rates of folding and unfolding should follow the Arrhenius equation over a limited range of temperatures. This study shows that molecular dynamic simulations of the unfolding of crambin between 500K and 560K do follow the Arrhenius equation. They also show that while there is a large amount of variation between the simulations the average values for the rate show a very high degree of correlation.


Introduction
Molecular dynamics (MD) simulations have become an important tool in understanding chemical and biochemical processes at the molecular level. Through the use of Newtonian mechanics and an empirically derived force-field, simulations have been used to investigate the interactions of drug molecules and their targets as well as to predict the behaviour of proteins and peptides (Vasquez et al., 1994). Originally simulations were carried out in vacuo, but now with the increasing power of computers, simulations are usually carried out using periodic boundary conditions in water.
Increasingly realistic simulations lead to a better understanding of processes such as protein folding and unfolding at the molecular level (Lindorff-Larsen et al., 2012;Scheraga et al., 2007). The literature on molecular dynamics simulations for protein folding and unfolding is extensive. This can include very large simulations that can be simplified using coarse-graining, down to simulations of small proteins or peptides. One area where MD simulations have been particularly widely used is in the simulation of prion proteins and the protein mis-folding diseases associated with them (Shamsir & Dalby, 2005;Shamsir & Dalby, 2007;van der Kamp & Daggett, 2010;van der Kamp & Daggett, 2011).
Molecular dynamics simulations give us a detailed view of protein folding and unfolding pathways and their rates of unfolding (Daggett, 2002). Temperature is often used to accelerate protein unfolding and it has been shown that this does not affect the protein unfolding pathway (Day et al., 2002). It is often difficult to relate the simulated results to experimental data. A review of simulated folding times has shown that the times predicted by MD and the experimental rates for thermal unfolding are in good agreement (Snow et al., 2005). Atomic force microscopy data is another possibility and this has been used to investigate protein folding of T4 lysozyme (Peng & Li, 2008).
The rate of protein folding at increasing temperatures should be described by the Arrhenius equation over a limited range of temperatures (Alberty): Where k is the rate of the reaction, A is a constant (pre-exponential factor) E a is the energy of activation, T is the Temperature in Kelvin and R is the Universal Gas Constant.
A rearrangement of the Arrhenius equation taking natural logarithms gives the linear function: A plot of the natural logarithm of the rate ln(k) against 1/Temperature will be a straight line if the simulations obey the Arrhenius equation.
The Arrhenius equation has been used in solid-state chemistry calculations but currently no studies have tested whether it is valid in MD simulations of protein folding (Huwe et al., 1999). This paper presents a MD simulation study using a small protein (crambin) to test whether the models do agree with the predicted linear behaviour.

Materials and methods
Crambin was chosen as the model protein for simulation of its size as it only has 46 amino acids (Caves et al., 1998). A high resolution crystal structure of crambin (3NIR.pdb) was downloaded from the RCSB Protein Databank (Rose et al., 2011;Schmidt et al., 2011). Molecular dynamics simulations were carried out using Gromacs 4.6.4 on an Ubuntu 12.04 machine with GPU acceleration (Hess et al., 2008).
The protein was solvated using periodic boundary conditions and a surrounding distance of 4nm around the protein. Simulations were run at 500K to 560K at 10K intervals using the OPLS forcefield. At temperatures of 570K and above the simulations fail to complete. The models were initially equilibrated using the canonical (nvt) and isothermal-isobaric (npt) ensembles. At 500K the simulations were run for 10ns with a time-step of 2fs. At 560K the simulations were run for 1ns with a time-step of 2fs. Secondary structure was calculated using DSSP (Dodge et al., 1998) (the data is available in dssp_files.zip). RMSD deviations from the original crystal structure were calculated in Gromacs and displayed in Grace (the data is available in rmsd_grace_datafiles.zip). All of the scripts for equilibration dynamics and analysis can be found in the accompanying data files (the Gromacs files are in gromacs_files.zip and the shell script to run all the simulations and analysis is gromacs_runs_ complete_analysis.sh).
The statistical analysis of the rate data was carried out using SPSS version 22 (IBM_Corp., 2013). The line of best fit to the mean data was fitted using linear regression (the data files are available as md_arrhenius_crambin_averaged.sav, md_arrhenius_crambin_ averaged.spv). The line of best fit for the complete dataset was calculated using the general linear model (the data files are available as md_arrhenius_crambin.sav, md_arrhenius_crambin.spv). All of the scripts and data for the statistical analysis are available in the supplementary materials.

Results
Over these time periods crambin did not unfold as much as had been expected. As it is such a small protein it seems to be particularly stable to rises in temperature. There were three possible end points that could be used in determining the rate: 1) The unfolding of the C-terminal final bend ( Figure 1, the green region from residues 36-38).
3) The increase of the root mean squared deviation (RMSD) to 0.4nm from the crystal structure.
It was not clear which of the measures would be the most reliable and so time was taken for all three. There were however missing values because the end points were not reached during the simulations. This was particularly apparent in the 540K simulations, which seem to be anomalous as can be seen in the boxplots for the reaction times from the simulations (Figure 2A-C). The boxplots show the expected downward trend, although there is considerable variation between the times taken for the repeated simulations. This variability declines at higher temperatures as rates become faster. A summary of the times to the three different end points are given in Table 1-Table 3.
The Arrhenius plot can be constructed using either the mean values for the different temperatures or all of the values for the repeats. In both cases a linear model is produced but the correlation is much stronger for the averaged values, which also give a narrower confidence interval for the line of best fit ( Figure 3A-C). The wider confidence intervals for all of the data and the extent of the scatter can be seen in Figure 4A-C.   The parameters for the lines of best fit using all of the data and only the mean values are given in Table 4 and Table 5.

Discussion
All three of the end points used produce a linear model in the Arrhenius plot with a high degree of correlation. This suggests that in these cases the MD simulations are following Arrhenius behaviour and that these models can be used to predict rates of unfolding.
Of the three end points the RMSD variation is the easiest to calculate and least ambiguous, but it is also difficult to understand what this signifies at the protein level, when compared to using the disruption of secondary structure as a metric. The other advantage of using the unfolding of the secondary structure is that experimental values are available for unfolding of proteins and so if an appropriate end point can be found in the simulations then the gradients of the Arrhenius plot can be used to calculate the activation energy for unfolding.
This study also highlights the high degree of variability between the trajectories of different simulations. This variability is very high at lower temperatures. Simulations were repeated ten times and this resulted in values for the standard errors for the time of unfolding that could be up to 25% of the time. These standard errors are large A C B Figure 4. A: The Arrhenius plot for the unfolding of the final bend using the complete data. B: The Arrhenius plot for the unfolding of the beta sheet using the complete data. C: The Arrhenius plot for the increase of the RMSD by 0.4nm from the crystal structure using the complete data.  and so the number of simulations that are run needs to be increased in order to reduce them. As the standard errors falls with the square root of the sample size this would mean a 4-fold increase in the number of simulations is needed to reduce the standard error by half. This would suggest that more reliable results could be obtained by carrying out 40 simulations at each of the temperatures, which is a considerable additional computational burden.
The large amount of variation also affected the quality of the lines of best fit to the Arrhenius equation. The coefficient of determination was much lower when considering all of the data but nonetheless there is clear evidence of the simulations following Arrhenius behaviour. The confidence intervals for the linear models using all the data and the average data are similar. The variation in the gradient is too large for making comparisons with experimentally derived energies of unfolding and this is another reason why a larger number of simulations will be needed in future studies.
There was a surprising degree of agreement in the slopes and intercepts of the bend loss and beta sheet loss end-points. This suggests that the energies involved in stabilising the beta sheet and final bend are similar. This consistency is encouraging and suggests that detailed energy predictions will be possible from MD simulations.
Crambin was not an ideal case for using in this study. Although it is very small and allows the simulations to be run in a shorter time, the protein does not unfold very much and so longer time-scales are needed within the simulations. Prion protein is another possible model that could be used to test unfolding as this has been shown to unfold over short simulation times (< 10ns) (Shamsir & Dalby, 2005). The other alternative is longer simulations times in order to produce clearer and less ambiguous end-points for the simulations. Author contributions ARD conceived the study and designed the experiments. ARD and MSS carried out the simulations and the analysis. ARD and MSS were involved in all of the versions of the manuscript and have agreed to the final content.

Competing interests
No competing interests were disclosed.

Grant information
The author(s) declared that no grants were involved in supporting the work.  . 1997 ) and as the current simulations are run at high temperatures where the Sep 30; 94(20): 10636-10640 protein is unstable, this is an issue. The dielectric behaviour of the water model at these high temperatures is also a concern for me.
I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
No competing interests were disclosed.

Dmitry Nerukh
Institute of Systems Analytics, Aston University, Birmingham, UK I have two main critical points: The conceptual one is that the loss of the beta sheet or the bend could not be, strictly speaking, classified as "unfolding". There is very little in the text that describes how specifically the authors quantified the moment of this loss of structure. Looking at Fig. 1 it can be seen that both motifs come back for some periods of time. Is this a loss of structure or just equilibrium fluctuations? If the latter, then there are no observable unfolding event in the data.
As authors rightfully admit, the uncertainty of these loss of structure times is very high, especially at 2. As authors rightfully admit, the uncertainty of these loss of structure times is very high, especially at low temperatures. This is most likely because the length of simulations, 1-10ns, is too short (as again, the authors mention this in the text) and also because the algorithm of determining these times is not statistically robust. For proteins, even small ones like crambin, typical time scales for structural rearrangements is of the order of tens of nanoseconds. Even dialanin takes hundreds of nanoseconds to accumulate statistically sound number of its very simple conformational transitions.
Taking into account these points, the main conclusion of the authors based on data shown on Fig. 2-4 appears unsubstantiated to me.
Overall, the paper is interesting and worth indexing if the authors more rigorously calculate the "loss of structure" times. I would consider building a statistical network on the major conformational states, in the spirit of Markov States Model. Then, the average transition times would faithfully represent the "unfolding" events.
I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
No competing interests were disclosed. Competing Interests: In the article, there is a comprehensive explanation of study design, methods and analysis. The conclusions are well explained and justified on the basis of the results, Furthermore, the authors have made the data fully available, including the scripts to reproduce their work, being able to be revised by the author. That point makes the work more trustable and robust.
As a consequence of that, the manuscript is recommended to be approved. However there are some minor changes and comments that the authors may consider to improving the study.
The authors stated that "The Arrhenius equation has been used in solid-state chemistry calculations but currently no studies have tested whether it is valid in MD simulations of protein folding". This phrase should be rewritten and clarified, because as it is written, it seems that this is the first time that the Arrhenius equation is used in protein folding simulations and that's not true because its usage is widely established in protein folding.
To gain statistical significance, the authors performed 10 replicas, which is a good number. However more replicas changing the force-field, for instance using Amberff99SB-ILDN, would be a good option to gain more statistics and at the same time ensure to avoid the force-field dependence of the results, as in folding/unfolding processes it is very relevant (see for instance ). Simulations Lindorff-Larsen ., 2012 et al with one or two more force-fields, performing 10 replicas for each one, would be appreciated.
To explore if the protein is unfolded or not, three different metrics were employed. One of them is the increase of the root mean squared deviation (RMSD) to 0.4nm from the crystal structure. Why 0.4 and not another value? There are some references to that in the literature or it is just an observed trend in the simulations? If it is an observed trend, would change varying the force-field. In , they Shimada , 2001 et al. identified important structural regions of Crambin. Some of them are used by the authors as unfolding indicators, although the aminoacid numbers are not exactly the same. Thus, I wondering if helix 1 (residues 6-18,sequence: SIVARSNFNVCRL) could be also a good indicator of Crambin unfolding with the employed PDB structure, and if it is, if there is a reason to not be used.
As the authors stated longer simulations of 40ns or even longer should show a more realistic and picture of Crambin unfolding. It is true that it constitutes a considerable additional computational time, but will probably result in stronger and more sound results. Thus I encourage the authors to do it, at least, in the future.

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
No competing interests were disclosed.

Competing Interests:
Author Response 30 Sep 2015 , University of Westminster Andrew Dalby I would like to thank the referee for his helpful comments.
We definitely need to correct that omission in referring to other work using Arrhenius on proteins, although it has mostly been on short segments of protein or peptides rather than a complete protein. In writing the introduction I was surprised by the lack of literature on using Arrhenius in a biological simulation context, except for enzyme simulations. I had searched for Arrhenius behaviour in molecular dynamics and as stated I only found references to its use in simulations of materials. It seems that the problem was I should have searched just for Arrhenius and in using behaviour it needs to be the US spelling, behavior to return the protein modelling literature! The work of Baker, Karplus and Kuriyan, and Pande are particularly important in establishing the field. The more recent work by Best and Mittal is another excellent example which incorporates the effects of forcefields.
I agree with the point about the forcefield but this adds a confounding variables and also increases the variance of the simulations, which increases the uncertainty in the gradient and has a negative impact on the confidence intervals for predicted energy barriers. An alternative is to use the technique to calculate relative differences in energy barriers between protein variants using the same forcefield. As you are calculating a difference the forcefield effects can be considered to cancel out. This was the idea behind free energy perturbation calculations.
I am the author Competing Interests: