Model

Abstract

Molecular dynamics simulations were performed in order to find out the heat stability of our fusion constructs. The simulated parts of our constructs were the ones for which we did not have any data from literature regarding the heat resistance. In particular, our red fluorescent protein and keratin binding domain were simulated between 298 K and 358 K. The goal of the simulations was to help us choose the best purification strategy once we were able to produce the fusion proteins. Based on the results, we decided to use a purification protocol where the proteins are heated up to 343 K, which is faster and cheaper than using His-tags and Ni-columns for purification.

Introduction

While we were designing our silk repeat containing fusion proteins, we had been provided with a purification protocol for a similar silk based fusion construct by one of our advisors. This protocol requires heating up the mixture of unpurified proteins to denature most of the non-silk based proteins. Therefore, it was necessary to see whether or not our modified fusion proteins would be able to survive this heat treatment. If they remained functional, we could then use this faster and cheaper protocol for purification compared to using His-tags and Ni-column chromatography. Also, due to the high complexity (e.g. a high GC content) of our silk repeat sequences, we had to wait approximately one month to get the synthesized parts. To save time, instead of experimental tests we decided to use molecular dynamics (MD) simulations to find out the heat resistance of our red fluorescent protein (mRFP) and keratin binding domain (KBD), for which we were unable to find the relevant data in literature. Our other parts were known to be heat stable so they were not simulated. For example the cellulose binding domain was already known to be heat resistant because it was from a thermophilic bacterium, Clostridium thermocellum.

Instead of separately simulating all of our fusion construct combinations that contained silk with either mRFP, KBD, or both, it was assumed that due to the linker sequences the behaviour of the parts that were simulated would be relatively similar to how they would act as individual non-fusion proteins. This assumption also saved a lot of computational resources and therefore time that would have been needed in order to simulate the much larger and higher amount of complete fusion proteins.

Methods

MD simulations were performed using GROMACS [1] software version 2016.4. The results were visualized using VMD [2]. The computational resources were from CSC IT center for science, Finland.

Prior to running the simulations, the structural files for both mRFP (PDB ID: 2VAD) and KBD (PDB ID: 1LM5) caused some adjustments to be made in order to successfully simulate these molecules. For mRFP, the main issue was that the chromophore, i.e. the part of the molecule that causes its color, is formed by three neighboring amino acids in the peptide chain reacting with each other during the maturation of the protein. The final protein structure therefore contains a non-standard amino acid, which means that the force field files need to be modified to contain the correct parameters for the chromophore. Obtaining and validating the parameters can be difficult, but luckily Dmitrienko et al. [3] had already calculated these for the OPLS-AA force field. This is also why all of our simulations were performed using the OPLS-AA force field [4] after manually adding the correct parameters for the chromophore. Some internal residues were missing from the KBD-C structural file, which is why the SWISS-MODEL fully automated protein structure homology-modelling server [5] was used to approximate the coordinates for the missing residues.

After cleaning up the structure files and making the necessary changes described above, the simulations were set up by creating a periodic box around the molecules so that there was at least 1.0 nm between the protein and the box sides. Three replicates were made for both mRFP and KBD-C. This was followed by filling the simulation box with water molecules using the TIP4P water model [6]. Na+ and Cl- ions were added to neutralize the charge of the protein. Energy minimization was performed, followed by a 0.5 ns NVT equilibration and a 10 ns NPT equilibration. The final production run contained 5x 45 ns time intervals at increasing temperatures (298 K, 313 K, 328 K, 343 K, and 358 K to be precise). The temperature was kept constant for the entire time intervals except for the last 100 ps of each interval during which the system was heated up by 15 K. 343 K was the main temperature of interest due to the purification protocol that was going to be used, but the other temperatures were simulated to see if there is any change in behaviour while the system is heated up and to make the rise in temperature more even.

Results

The aim of the simulations was to find out if our proteins would denature when heated up to 343 K. Therefore, the analysis was focused on seeing if there is any change in the structure of the proteins that could indicate potential denaturing. Thus, the radius of gyration of the proteins as well as the number of hydrogen bonds within the proteins were calculated at each time interval. The first 5 ns of each time interval was not analysed to reduce possible errors caused by the the system getting used to the sudden increase in temperature. The data points and error bars were calculated by averaging the results from each time point within the time intervals and from each of the three replicate simulations.

The radius of gyration (Rg) describes the dimensions of a protein by calculating the mass weighted distances of the parts of the protein around its center of mass. Significant changes of Rg at higher temperatures could mean that the tertiary structure has changed enough to destroy its functionality. The average radii of gyration as a function of temperature for mRFP and KBD-C are shown in Figure 1.

Figure 1. Radius of gyration as a function of temperature. Red circles and error bars are for mRFP, black triangles and error bars are for KBD.

As can be seen from Figure 1, the radius of gyration does not change significantly at any temperature for either protein. This may indicate that the proteins can tolerate these temperatures without denaturing. To further test this, the number of hydrogen bonds within the simulated proteins was analysed. Hydrogen bonds play a major part in the formation of both the secondary and tertiary structures of proteins. Therefore, a significant decrease in the amount of hydrogen bonding between amino acids would indicate denaturing of the proteins. The number of hydrogen bonds is shown as a function of temperature in Figure 2 for mRFP and KBD.

###picture here###

Figure 2. The number of hydrogen bonds between amino acids in the simulated proteins as a function of temperature. Red circles and error bars are for mRFP, black triangles and error bars are for KBD.

As can be seen from Figure 2, there is a small decrease in the number of hydrogen bonds in mRFP, whereas there is no noticeable difference in KBD. A small decrease can be expected, because higher temperatures translate to more movement of atoms.

In conclusion, based on the results both mRFP and KBD-C should be able to stay intact at 343 K. This allowed us to purify our silk containing fusion proteins with the heat treatment protocol provided to us instead of using His-tags and Ni-column chromatography, saving time and resources.

References:

1 Berendsen, H. J. C., van der Spoel, D., van Drunen, R. GROMACS: A message-passing parallel molecular dynamics implementation. Comp. Phys. Comm. 1995, 91 43–56.

2 Humphrey, W., Dalke, A., Schulten, K. VMD: visual molecular dynamics. J Mol Graph 1996, 14(1), 33–38.

3 Dmitrienko, D.V., Vrzheshch, E.P., Drutsa, V.L., Vrzheshch, P.V. Red fluorescent protein DsRed: Parametrization of its chromophore as an amino acid residue for computer modeling in the OPLS-AA force field. Biochemistry (Moscow) 2006, 71(10) 1133-1152.

4 Jorgensen, W. L.; Tirado-Rives, J. The OPLS Potential Functions for Proteins. Energy Minimizations for Crystals of Cyclic Peptides and Crambin. J. Am. Chem. Soc. 1988, 110 (6), 1657–1666.

5 Waterhouse, A., Bertoni, M., Bienert, S., Studer, G., Tauriello, G., Gumienny, R., Heer, F.T., de Beer, T.A.P., Rempfer, C., Bordoli, L., Lepore, R., Schwede, T. SWISS-MODEL: homology modelling of protein structures and complexes. Nucleic Acids Res. 2018, 46(W1), W296-W303.

6 Jorgensen, W. L.; Madura, J. D. Temperature and Size Dependence for Monte Carlo Simulations of TIP4P Water. Mol. Phys. 1985, 56 (6), 1381–1392.