Asymmetric Behavior of Thymidylate Synthase Dimer Subunits in Denaturating Solvent Observed with Molecular Dynamics

A molecular dynamics simulations of the thymidylate synthase denaturation in chaotrope solvents (urea, guanidinium hydrochloride) were performed on 600 ns timescale. It appeared that this dimeric enzyme undergoes partial unfolding asymmetrically. It was shown also that urea is a better denaturant in the MD condition, as compared to guanidinium chloride. The unfolding occurs first at the external helices (AA 88-118) and follows by the AA 188-200 region. The present results correspond to the suggested in the literature activity of thymidylate synthase through a half-the-site mechanism.


Introduction
The thymidylate synthase enzyme has been studied experimentally for a long time due to its unique role as an agent catalyzing reductive methylation of dUMP using 5,10methylenetetrahydrofolate as a cofactor [1].A distinctive feature of thymidylate synthase is that it is the sole de novo source of thymidylate in a cell, which makes it the target protein in various types of chemotherapies [2][3][4].
Thymidylate synthase (TS) is a dimeric enzyme.There is also an interesting feature of TS.Two structurally equivalent active sites of TS display apparent functional nonequivalence suggesting TS to be a half-the-sites activity enzyme [5].The principle of half-the-site reactivity was established by Levitzki and coworkers [6] summarizing the discoveries in the late 1960s of many enzymes exhibiting this phenomenon.By this term they meant that the proteins composed of identical subunits with  potential active sites, reacted with a substrate or an inhibitor so that only /2 sites were occupied when enzyme was saturated with the ligand.
TS undergoes unfolding and denaturation under the action of chaotropes like urea or guanidine hydrochloride (also guanidinium chloride, usually abbreviated as GdmCl or GuHCl).The early studies have shown that chaotropeinduced unfolding can result in the formation of monomeric species under the action of about 6 M urea and 4 M guanidine hydrochloride solutions [7][8][9].Asymmetric behavior, if observed under denaturating conditions, could improve understanding of the half-the-site mechanism.
Chaotrope denaturation mechanism was thoroughly examined for polipeptides and proteins using a computer simulation method, the molecular dynamics (MD) [10][11][12][13][14][15][16][17][18][19][20].In this method time evolution of an atomistic structure of a protein and surrounding solvent molecules is studied, by solving Newton equation of motions.Potential energy in MD is defined by a so called force field, which contains equations and parameters to find the energy.Applicability range of MD is however limited by currently available hardware resources.
Aim of the already performed chaotrope MD denaturation studies was describing transition states and unfolded states.Achieving it with standard MD protocols was too expensive, so a way to accelerate the method had to be found.Some of the referenced studies use elevated temperatures or increased denaturant concentrations to allow system to quicker crossing of conformational barriers and in result faster unfolding.There are also studies using enhanced sampling [20], where biomolecular system is simulated in parallel in multiple temperatures or states (different level of unfolding).These different states can be exchange under certain conditions, allowing much better probing conformational space exploration.
These methods come however with a price of lowering simulation fidelity.Our work is focused on early steps of unfolding and we wanted to describe sequence of events that starts denaturation, rather than getting into unfolded state.Therefore we have decided to use a standard MD protocol, without unphysiological temperatures or enhanced sampling.
In this work we present a study of TS enzyme using MD in water and chaotrope denaturants.We show that denaturation conditions result in asymmetric unfolding of the units and on urea example we give insight into mechanism of first step of TS unfolding.

System Preparation.
Thymidylate synthase from rat crystallographic structure was taken from the PDB database (PDB ID: 2TSR [1,21]).The structure was resolved using X-ray diffraction at 2.6 Å resolution.Crystallographic structure is a dimer of dimers, that is, it consists of two equivalent copies of the enzyme.These copies are later referred to as AB (first dimer, chains A and B) and CD (second dimer, chains C and D).The structure further used in MD simulations was stripped of ligands (Tomudex drug and dUMP), however crystal water molecules were preserved.
Solvent mixture environments for the MD simulations were prepared using a Packmol application [22].Boxes of 80 Å × 80 Å × 80 Å in size were filled with solvent molecules according to expected concentrations.Required densities for 8 M urea and 6 M guanidinium chloride solutions in 310 K (37 ∘ C) were calculated using the empirical equations presented by Kawahara and Tanford [23] as well as tables with water densities depending on temperature.
For a model of 8 M urea solution, the box was composed of 10963 water molecules and 2467 urea molecules; the values proceed from calculated density of 1.1132 g/cm 3 and correspond to 0.750 grams of urea per gram of water, which matches well to the experimental value (0.755) [24].With calculated GuHCl solution density equals 1.1363 g/cm 3 , a box containing 1850 guanidinium ion molecules (Gu + , accompanied by chloride anions) and 9638 water molecules was generated, corresponding to 1.018 grams of denaturant per gram of water, which is in good agreement with the experimental value (1.009) [24].The boxes were subsequently minimized in NAMD [25] for 500 ps, followed by heating (from 0 K to 310 K), and proceeded by 10 ns molecular dynamics simulation.
The PDB model of dimeric protein of rat thymidylate synthase was immersed in each of the prepared boxes and all solvent molecules within 1.5 Å of the protein were removed, excluding crystal water.As water is much smaller in size than guanidinium ion/urea molecule, the aforementioned procedure of solvent removal was biased towards lowering non-aqueous solvents concentrations.The system was later neutralized with sodium cations, since TS in these pH conditions is negatively charged; also removal of some Gu + ions during protein insertion needs to be complemented by cations.Six simulation systems were created: simulation of AB dimer in water, urea, and guanidinium chloride, as well as CD dimer in water, urea, and guanidinium chloride.We have performed independent simulation with two different dimer structures to reduce bias from starting coordinates.

Molecular Dynamics Setup.
During the MD simulation the system was described with the use of the CHARMM22 force field [26,27], water molecules were represented by the three-site TIP3P model [28], urea parameters were characterized by Caflisch and Karplus [29], and guanidinium ion parameters were taken from Kuczera [26]; parameters for both nonaqueous solvents were consistent with the CHARMM22 force field.Periodic boundary conditions were used with the Particle Mesh Ewald method [30] to account for long-range electrostatic interactions, with grid spacing of less than 1 Å.Simulations were performed using the NAMD 2.8 [25] engine on an IBM BlueGene/P supercomputer.
The minimization of the system internal energy was performed in three steps, each for 100,000 steps: (i) only solvent and ions, (ii) solvent, ions and C  protein atoms, and (iii) the full system.This was followed by actual MD simulation with 1 fs time step: heating, equilibration, and production.Heating was performed in constant volume conditions until 310 K temperature was reached, 0.8 ps for every 10 K step, and giving 24.8 ps run and followed by a 3.2 ps equilibration run at 310 K.These steps were followed by a production run of hundreds of nanoseconds, performed in NPT ensemble at 310 K, with temperature stabilized using a Langevin thermostat (dumping coefficient of 5 ps −1 ) and pressure set using a Langevin piston [31,32].Langevin method used for keeping temperature and pressure constant are stochastic methods, introducing random fluctuations to the system.Therefore subsequent trajectories of the same system might differ [33].Accounting for this fluctuation is possible with having all simulations in all the solvent conditions performed with either AB or CD dimer.

Denaturant Concentrations.
Setting pressure constant in the simulation accounts for imperfections of structure preparation.If any cavities were left after inserting the protein to a solvent box, the periodic boundary conditions box will be compacted to account for that.However, due to variable volume, the simulation is in principle not a constant concentration simulation (number of solvent molecules do not change with volume fluctuations).
Taking this into account, as well as effects of insertion of the protein to the box effective concentrations in this study, come to ca. 5.9 M for urea and 3.8 M for GuHCl, at the beginning of simulation.Being consistent with experimental values of 6 M urea and 4 M GuHCl provided in the Introduction.

Trajectory Analysis.
The trajectories obtained during the simulation were analyzed using the Cpptraj application, available in the AmberTools 13 suite [34,35].The following aspects of the simulation were taken into consideration.
RMSD value is minimized over all possible translations and rotations of the superimposed molecule.In this paper we consider RMSD calculated only for C  atom positions.
(ii) Intermonomer RMSD: A measure of asymmetry between dimer units of the protein in function of time.In each time step RMSD value is calculated between the first and the second monomer unit coordinates (coordinates of the second monomer at time  substitute reference position in (1)).Coordinates are translated and rotated using Kabsch algorithm [36] before RMSD calculation to minimize RMSD value.
(iii) Root Mean Square Fluctuation (RMSF): A measure of mobility of atoms/residues during a period of time.
If x  () is the th atom position (out of ) after  time steps (out of ), and x ref  is the th atom position reference value, then RMSF is defined as As in case of RMSD, in this paper, we report RMSF of only C  atoms to remove contributions from internal movements in amino acids.The reference value of atom position is its average position over time.Structures were superimposed to the reference (starting) structure according to the lowest RMSD before calculating RMSF.
(iv) Solvation Number: The number of solvent molecules occupying space in the 5 Å cut-off distance from the solute.

Effect of Solvent Choice on the Thymidylate Synthase Struc-
ture.Initially we have simulated all 6 solute/solvent pairs on 300 ns time scale.General observation from RMSD measurements (Figure 1) is that the protein in water remains stable, so it can be taken as a reference for chaotrope simulations.There are ongoing changes in conformation in the urea solution.
Guanidinium hydrochloride solutions shows also signs of denaturation, but much slower.Progress of denaturation is slow, when compared with smaller molecule studies.
In terms of asymmetric behavior in Figure 1 we present intersubunit RMSD data for TS simulations.In water simulations the intersubunit RMSD remains at level of 2 Å.For guandinium hydrochloride in one simulation the value raised to more than 3 Å, but in other simulation remained at 2 Å level.Finally in case of urea the value raised to 5-6 Å level.

In-Depth Analysis of Thymidylate Synthase AB Dimer
Denaturation in Urea Solution.As 300 ns simulation of thymidylate synthase in urea revealed signs of ongoing denaturation, we have chosen to continue a longer simulation and we have been able to run MD for twice as long (600 ns) for the following systems: (a) AB dimer in urea, (b) CD dimer in urea, and (c) AB dimer in water.The trajectories allowed to perform an in-depth analysis of the obtained results to answer questions about initial steps of denaturation.
Comparison of the 600 ns trajectories is presented in Figure 2. Global trend for the three trajectories is not different from that of the 300 ns simulations.The protein remains stable in water while denaturation in urea continues, as seen on the RMSD plot.
Also intermonomer RMSD continues to rise to level of approximately 10 Å by the end of a simulation time.This value is higher than RMSD of the whole protein calculated from starting structure, in range of 6 Å by the end of simulation.On the other hand there is no visible sign of increasing solvation of the protein, as seen on model systems simulations [20].
To approach the problem profoundly, local measures have to be taken into account.Figure 3    Since the simulation in water did not show signs of structure evolution, the whole trajectory of the protein in water is taken as a reference.Residues are numbered according to the protein sequence, with certain residues not present in the PDB structure (i.e., amino acids 1-20 and 302-307).
during the simulation of the protein in urea, as compared to simulation in water.Figure 4, on the other hand, presents how secondary structure of the protein is affected by denaturing conditions.The regions with the highest observed flexibility are -helices between amino acids 88 and 138.In both A and B monomers one can observe decrease in percentage of time in which amino acids are in helical conformation (see Figure 4).General picture is similar, however careful observation of atoms mobility reveals that between 400 and 500 ns subunit A shows approximately 1.5 higher values of fluctuations than the B one (see Figure 3: level of 12-13 Å for the monomer A and maximum value of 8 Å for the monomer B).On the other hand, the region between amino acids 138 and 150 is occupied in the monomer A by anti-parallel strands, which are not present in the monomer B. Careful analysis of Figure 4 might even suggest that these anti-parallel strands in the B subunit are slowly developing in the progress of the simulation.In case of the A subunit there is also a small increase in fluctuations in Figure 4: The evolution of the secondary structure of the AB dimer of thymidylate synthase in urea solution calculated using DSSP [38] algorithm for different intervals of the trajectory.Presented are helical areas (red) and  strands (blue).Darker color represents secondary structure element present in at least 75% of frames in the interval; lighter color represents element present in 25-75% of frames.Residues are numbered according to the protein sequence, with certain residues not present in the PDB structure (i.e., amino acids 1-20 and 302-307).
of the most flexible regions in the course of initial steps of denaturation are visualized in Figure 5.

Conclusions
In this work we have shown that denaturation of a relatively large protein, like thymidylate synthase, in adverse solvent conditions, might be observed in molecular dynamics simulations on hundreds of nanosecond timescales (0.5 s).Shorter simulations, requiring less computational effort would not be able to catch even a limited scope of the denaturation.The fact that the protein remained stable in similar simulation setup, but only with pure water, allows to show that urea and guanidine hydrochloride used in simulations are responsible for denaturation effects seen in Figures 1-4.However early signs of denaturation where seen in the simulation, simulation time shorter than a microsecond is not enough to capture full denaturation pathway of the protein.Similar timescales were successful to observe denaturation for small polipeptides, like a lysozyme [19], which has more than 5 times less residues than the thymidylate synthase dimer.Even case of such smaller proteins only partial denaturation pathway was observed during 1 microseond.Therefore, as expected our trajectory does not give answers on unfolding pathway of the internal core, including the active site.
When analyzing denaturation of the protein monomers only denaturation of external loops can be seen in the trajectories.According to the results, unfolding of the external helices comprising the amino acids from 88 to 118 is the first step in the denaturation process of thymidylate synthase, later accompanied by unfolding of the 188-200 region.Since the latter region is close to the active site, we suspect that this is where denaturation of the active site is beginning.Preferential location of 88-118 amino acids as well as 138-150 ones on the exterior of the protein (see Figure 5) makes them ideal areas for denaturation to start at.There is no surprise they show very high atoms mobility (see Figure 3), followed by losing their secondary structures in the progress of the simulation (see Figure 4).
The denaturation process ran asymmetrically with the greater changes and faster unfolding appearing in one subunit rather than in its counterpart.This is similar to the proposed mechanism of the dimer activity, where only one monomer is active at the time [37].The results we obtained do not prove the half-of-the-sites activity of TS, rather they show unequivalent unfolding of monomeric subunits which may form a molecular basis for the half-reduced activity of the dimer.

Figure 1 :
Figure 1: Root mean square deviation in reference to starting structure (black) calculated between C  carbon coordinates of the first dimeric subunit versus the second subunit (intermonomer RMSD red) of the thymidylate synthase in function of time for 300 ns molecular dynamics simulation.Raise of the intermonomer RMSD value shows asymmetric unfolding of the subunits, especially in case of urea.
(i) Root Mean Square Deviation (RMSD): A measure of similarity of molecule 3D-fold to a reference structure, usually crystallographic one, is routinely used in analyzing MD simulations to find stability of a molecule in time.If x  () is the th (out of ) atom position after time , and x ref  is value of the atom position in a reference molecule, then RMSD is defined as

Figure 2 :
Figure 2: Comparison of RMSD (left side, black), intermonomer RMSD (left side, red), solvation number of water (right side, black), and urea (right side, red).Values were processed with the moving average method (window of 20 points) to decrease noise.

Figure 3 :
Figure3: RMSF of the protein in urea solution during different parts of trajectory.Black line corresponds to RMSF of the first monomer (A or C) and red line corresponds to RMSF of the second monomer (B or D).Since the simulation in water did not show signs of structure evolution, the whole trajectory of the protein in water is taken as a reference.Residues are numbered according to the protein sequence, with certain residues not present in the PDB structure (i.e., amino acids 1-20 and 302-307).

Figure 5 :
Figure 5: (a) Thymidylate synthase dimer 3D structure, colors represent 2D folds: red--helix, yellow--sheet, green-loops, blue molecule-ligands in the active site (not present in MD), and red beads-water.(b) TS structure with residues 88-118 (in both dimeric units) highlighted in black.(c) TS structure with residues 138-150 (in both dimeric units) highlighted in black.(d) TS structure with residues 180-200 (in both dimeric units) highlighted in black.
presents RMSF changes 3 10 -helix between residues 180 and 200.This region is an inner part of the protein, close to the active site.It is possible that large fluctuations in the A88-138 region result in noticeable movements of A180-200 region, and since fluctuations in the B88-138 are smaller in magnitude, they do not affect B180-200 area.The positions