Conformational relaxations of human serum albumin studied by molecular dynamics simulations with pressure jumps

Aim. In this work we developed a novel technique of obtaining the spectrum of conformational relaxations in a solvated protein using non-equilibrium molecular dynamics simulation. Methods. Structural relaxations in the protein are initiated by the abrupt jump of pressure in the NPT simulations. The change of the protein volume after the jump is monitored and the Maximum Entropy Method is used for spectral decomposition of the volume relaxation curve. The human serum albumin (HSA) is used as a test case. Results. The obtained relaxation spectrum of HSA contains one component attributed to the bulk water and five components caused by the relaxations of the protein globule and its hydration shell. All relaxation components are in good agreement with the available experimental data obtained by the time-resolved spectroscopy and the broadband acoustic spectroscopy of HSA. Conclusions. The developed technique allows obtaining spectra of conformational relaxations of soluble proteins in a range of time scales from ~0.1 ps to ~50 ns utilizing single non-equilibrium molecular dynamics simulation.


Introduction.
It is well known that the dynamics of protein globules extents from femtoseconds to seconds and even minutes [1].It is not surprising that none of existing experimental techniques is able to study the whole spectrum of motions in proteins.Various time-resolved spectroscopic methods could detect motions at the picoseconds to nanoseconds time scale [2].Other techniques, such as broadband acoustic spectroscopy, can cover characteristic times from nanoseconds to microseconds and larger [3].The majority of available experimental techniques are based on measuring relaxations in protein globules initiated by various factors, such as excitation of the chromophore in optical spectroscopy or local changes of pressure in acoustic measurements.Although very different stimuli may trigger relaxation in the system the general picture of the protein dynamics remains mostly consistent in different experimental techniques.Spatial resolution of optical spectroscopy is limited by the position of chromophore, while acoustic spectroscopy is a completely integral method, which is unable to focus on particular parts of protein structure.Thus, relaxation components are quite hard to attribute to particular conformational changes in the protein globules.
Molecular Dynamics (MD) simulations are now con sidered as an important complementary method, which allows interpreting experimental data in terms of detailed atomistic model and mapping them into the protein structure [4].MD is particularly successful in reproducing the data of time-resolved spectroscopy by computing equilibrium fluctuations of atoms in the vicinity of chromophores [5][6][7].To our knowledge no attempt was made to obtain the general spectrum of structural relaxations of a whole solvated protein in MD simulations.Such spectrum would provide a global overview of protein dynamics at very different time scales and could be compared with the results of integral tech-niques, such as acoustic spectroscopy.The origin of relaxation components could be determined later and mapped to particular conformational changes in the globule.
The goal of this work is developing a technique for obtaining the spectrum of conformational relaxations of a whole protein globule in the wide range of time scales in MD simulations.Instead of computing the equilibrium fluctuations of atoms and deducing relaxation times from them we put the whole system into highly nonequilibrium conditions and record relaxations directly.
The Human Serum Albumin (HSA) was used as a test protein.Relaxation dynamics of this protein is studied by various experimental techniques such as broadband acoustic spectroscopy [8,9] and time-resolved optical spectroscopy [10][11][12][13].This makes it an ideal test object for our technique.
The major goal of this study is developing a consistent methodology of detecting volume relaxations in MD simulations rather then detailed study of a particular protein.
Materials and methods.Molecular dynamics simulations.S i m u l a t i o n s e t u p.All MD simulations were performed with the GROMACS suit of programs (versions 4.0.2 and 4.5.1)[14,15].The GRO-MOS force fields G53A6 was used [17].The long-range electrostatic interactions were computed using the fourth-order PME method with a Fourier spacing of 0.15 nm.The cut-offs of 1.4 nm and 1.0 nm were used for the short-range Coulomb and the Lennard-Jones interactions respectively.Bond lengths within the protein were constrained using the LINCS algorithm [18].The water molecules were constrained using SETTLE [19].The usage of heavy hydrogens [20] allowed the time step of 4 fs in all simulations.Initial structure of the human serum albumin was extracted from the PDB entry 1AO6:A.
The protein was solvated by 23504 SPC water molecules [21] in rectangular periodic box.Necessary number of the counter ions was added to neutralize the system.The protein and the water with counter ions were coupled to independent heat bathes at 300 K using v-rescale thermostat with a coupling constant of τ t = 0.1 ps.Berendsen barostat with a coupling constant of τ p =1.0 ps was used for isotropic pressure coupling.E q u i l i b r a t i o n.Initial equilibration of the system at pressure of 1 atm was monitored by the C α RMSD of the protein relative to the starting frame and the secondary structure content.Fig. 1 shows that equilibration takes about 100 ns before both parameters could be considered fluctuating around settled equilibrium values.
P r e s s u r e j u m p s e t u p.In order to initiate relaxations in the system the reference pressure of the barostat τ p was changed abruptly and the simulation was continued with the coordinates and velocities of the atoms, which were recorded before the jump.After that the simulation continues at new pressure and the volume of the simulation box relaxes to new equilibrium value.The relaxation of the volume of the protein globule dominates this process because the water is rather incompressible and relaxes quickly.Thus, analyzing the relaxation of the box volume one could extract information about the relaxation of the protein globule.The individual exponential components of this relaxation will correspond to particular structural rearrangements in the protein globule.This technique is related to the well known fluctuation-dissipation theorem [22].In order to make the relaxation detectable we increased the pressure from 1 atm to 10000 atm in MD.This does not lead to any significant changes in the secondary structure content, which convinces us that the protein is not unnaturally distorted by the pressure jump.
P r o d u c t i o n s i m u l a t i o n s.In order to resolve relaxation components, which are scattered in time over several orders of magnitude, the coordinates were saved with increasing intervals.For the interval from 0 to 1 ps the coordinates were saved each 0.002 ps.For the times from 1 ps to 10 ps each 0.02 ps.For the times from 10 ps to 100 ps 0.2 ps, etc.The longest part of trajectory from 100 ns to 175 ns was saved each 800 ps.Such setup provided enough data points for all covered time scales.The volume of the simulation box was computed for each saved frame.
P u r e w a t e r r e f e r e n c e.In order to distinguish the relaxation of the protein globule from the response of bulk water the reference simulations of pure SPC water were performed.The same number of water molecules presenting in the system with solvated protein (23504) was placed into the cubic box and simulated using the same protocol with the pressure jump for 10 ns.
The maximum entropy method.The major problem in the analysis of volume relaxations is decomposition of relaxation curves into exponential components.One of the most successful techniques of such decomposition is the Maximum Entropy Method (MEM) [23,24].MEM is widely used in such areas as time-resolved spectroscopy [24,25], neutron scattering [26] and image processing in astronomy [27] and tomography [28].In this work we provide the first application of MEM to the analysis of volume relaxations in proteins.
Let us consider the abstract problem of fitting decaying curve G(x) by exponential spectral components exp(-t/τ): where A(τ) is a continuous spectrum of experimental signal, which should be found, τ is the relaxation time.In practice both the signal and the spectrum are discrete, thus where N is the number of discrete spectral components used to approximate the real continuous spectrum.In order to find the weights of spectral components A j one should minimize the deviation between the real signal and its approximation given by (2) in respect to A j .This could be achieved by minimizing χ 2 function: ( ) where M is the number of experimental points in the signal.
It is well known that such multidimensional optimization problems are ill-posed [23,25].This means that there are many possible solutions, which provide the same accuracy in terms of the mean square deviation from given experimental curve.We follow simplified formulation of MEM of Steinbach et al. [25], which allows solving such ill-posed problem.An auxiliary functional Q is composed as where S(A) is the Shannon-Jaynes entropy, L is unknown Lagrange multiplier.The S(A) reads where A 0 is the initial guess for A. It is possible to prove that maximization of the functional (5) produces a unique solution of initial optimization problem [23,25].According to (4) S is maximized under the constraint imposed by the χ 2 of real signal G and its approxima tion G fit .Such maximization ensures that the resulting set of values A j possesses minimal information content among all sets, which are able to approximate the signal G with given accuracy.In the other words the solution is free from spurious correlations and the distribu tion of values in A j is as smooth as possible [23,25].
The following algorithm was used to maximize the functional (4): 1) uniform values are assigned to initial guess σ y 0 ; 2) initial value of L is chosen; 3) functional ( 4) is maximized using simple sequential unconditional optimization until the difference in the values of Q between two subsequent steps (tolerance) drops below given criterion; 4) L is increased; 5) go to step 3 until χ 2 drops below given value.Custom MEMfit software (http://members.multimania.co.uk.memfit/) was used to perform MEM computations.Detailed information about this software is available by request.Fitting of the relaxation curve of pure water was performed using N = 200, while fitting of the curve of solvated protein was performing with N = 400.
Results and discussion.Volume relaxation curves.The evolution of the periodic box volume after the pressure jump in MD simulations of solvated protein is shown in Fig. 2. The curve for pure water looks similar (data not shown).Apparently the curve shape is close to single exponential, which drops very fast in 0.01-0.1 ps.However, there are clear departures from a single exponential fit visualized in Fig. 2. It is clear, that this curve contains both very fast (~0.1 ps) and very slow (~50 ns) relaxation components.
Relaxations of pure water.Fig. 3 shows the spectrum of exponential components obtained by MEM from the volume relaxation curve of the pure water system.There are three well defined components visible as sharp peaks.The first peak τ 1 = 0.189 ± 0.12 ps is close to the relaxation time of the v-rescale thermostat used in our MD simulations (τ t = 0.1 ps), while the second peak τ 2 = = 0.592 ± 0.26 ps is close to the relaxation time of the Berendsen barostat (τ p = 1 ps).These two peaks are likely to be methodological artifacts.The third peak τ 3 = = 1.588 ± 0.33 ps is the only one, which could be attributed to the relaxation of water molecules.It is of the same order of magnitude as the time of dielectric relaxation of SPC water τ M , which varies from 3.5 to 4.1 ps depending on simulation parameters [29].
Relaxations of solvated HSA.Fig. 4 shows the spectrum of exponential components obtained from the volume relaxation curve of solvated HSA protein.This spectrum is much more complex then for pure water and contains eight peaks.The first three peaks τ 1 = 0.079 ± ± 0.046 ps, τ 2 = 0.521 ± 0.27 ps and τ 3 = 1.504 ± 0.66 ps are essentially the same as observed for pure water in Fig. 3. τ 1 and τ 2 originate from the thermostat and the barostat weak coupling algorithms respectively, while τ 3 is probably caused by the fast relaxation of bulk water.Other five peaks are absent in pure water and correspond to much slower relaxations of the protein and its hydration shell.These peaks are τ 4 = 7.53 ± 1.52 ps, τ 5 = 86.5 ± 8.4 ps, τ 6 = 545.9± 81 ps, τ 7 = 3289.4± 317 ps and τ 8 = 49787.6± 4803 ps.
Comparison with experimental data.We compared our data with the results of the broadband acoustic spectroscopy [8,9] and optical spectroscopy of the HSA [10][11][12][13].Four major components were detected in the acoustic spectra of HSA in various experimental conditions [9].The shortest relaxation time of 200-600 ps may reflect the motions of small flexible segments on the globule surface.Our relaxation time τ 6 falls into this region.The second experimental relaxation time is of order of 1-3 ns, which is close to our component τ 7 .The third experimental relaxation time is of order of 40 ns, which fits to our component τ 8 nicely and is usually attributed to cooperative segmental motions of the protein controlled by pH.The slowest experimental component lies in the microsecond region and can not be reached in our simulations.The data of the time-resolved fluorescent spectroscopy of HSA revealed several relaxation components, which are attributed to either motions of the protein in the vicinity of chromophore or the dynamics of its hydration shell.The phase-fluorimetry of dielectric relaxations in HSA revealed two relaxation components in the regions of 0.4-2.0ns and 4.9-8.1 ns [10], which correspond to our components τ 6 and τ 7 .Time-resolved spectroscopy of the single tryptophan residue W214 in HSA revealed relaxation times of 25-100 ps [13].Relaxations of the surface hydration shell assisted by fluctuations of the protein is estimated to occur at the time scale of 100-150 ps in several proteins including HSA [12].These relaxations correspond qualitatively to our component τ 5 .
The component τ 4 could be related to fast components of hydration dynamics, which occur at 3-5 ps time scale [12].Relaxation times of 5.6 ns and 41 ns were reported in other time-resolved fluorescence studies of W214 in HSA [11], which agrees reasonably well with the components τ 7 and τ 8 .
Studies of the protein relaxations in MD simulations are usually limited to computing equilibrium fluctuations.This technique is not suitable for determining the whole spectrum of relaxations, which occur in the protein globule after some global stress.In this case the direct observation of non-equilibrium relaxations is preferable because it provides the whole spectrum of relaxations present in the single MD trajectory.Such global overview of protein relaxations is useful for general understanding the protein dynamics and could be compared directly with experimental data.
The necessary prerequisite for the direct observation of relaxations in MD is robust and reliable technique of spectral decomposition of obtained decay curves.We utilized the Maximum Entropy Method, which selects the most plausible solution of decomposition problem with minimal information content.
Another prerequisite is a sufficient amount of data points for all time scales on the relaxation curves.We solved this problem by using variable sampling intervals increasing with an increase of time.
Finally, the last prerequisite is an adequate initial stress, which triggers relaxations in the system.The jump of pressure is a good candidate for such stimulus because it acts globally on all atoms in the system and initiates the whole range of relaxations with very different time scales.
Comparison of volume relaxations occurring in our MD simulations with the data of various experimental techniques is only possible on the semi-quantitative level.We managed to attribute six relaxation components to either experimentally observed relaxations in HSA or to the relaxation of bulk water.The nature of all observed relaxation components is summarized in Table .Components τ 3 and τ 4 are attributed to the motions with overlapping time scales, but there is no ambiguity because the component τ 4 is not observed in pure water and thus clearly corresponds to hydration dynamics of the protein.Despite quite uncertain semi-quantitative nature of the acoustic data for HSA [8,9] we managed to map three relaxation times τ 6 -τ 8 to experimental relaxation components in this technique.In general, it is remarkable that all relaxation components found in our simulations agree with the experimental data reasonably well.
To our knowledge the relaxation times of order ~50 ns were never detected before in MD studies.The sensitivity of our technique allows detecting even slower relaxations providing that MD trajectory is long enough.
Conclusions.The protocol of MD simulations, which allows obtaining the spectrum of conformational relaxations in the solvated protein at the time scales from ~0.1 ps to ~50 ns using a single simulation, is developed.The method utilized an abrupt pressure jump to create non-equilibrium stress in the system and uses MEM fitting of the volume relaxation curve to extract individual relaxation components.
The structural relaxations of HSA, obtained by different experimental techniques, are reproduced with good accuracy in our protocol.
Acknowledgments.The author is grateful to T. O. Hushcha for inspiring MD studies of HSA and discussing the results.V. N. Kharkyanen and N. M. Berezetskaya are acknowledged as developers and testers of the MEMfit software.

Fig. 1 .
Fig. 1.Evolution of the secondary structure content (A) and the C α RMSD with the starting structure (B) in the course of equilibration MD run

Fig. 2 .Fig. 3 .
Fig. 2. Evolution of the periodic box volume in the course of production MD run.Curve 1 shows single exponential fit.The pressure jump occurs at time t = 0. Insets show the regions of small and large times separately

Fig. 4 .
Fig. 4. The spectrum of relaxation times obtained by MEM from the volume relaxation cur ve of solvated HSA shown in Fig. 2. Positions of the peak maxima are shown.Log scale on X axis and different scale on Y axis are used after τ = 3 ps to show the peaks with large times and small weights