In this section, you will analyze two non-equilibrium properties of proteins, the diffusion of heat and a phenomenon known as temperature echoes. This section requires some knowledge of the principles of statistical physics, but the results of the simulations can be understood intuitively even without a deep understanding of their principles.
![\framebox[\textwidth]{
\begin{minipage}[r]{0.9\textwidth}
\noindent{\textbf{Ob...
...s of your simulation, the thermal diffusivity of your system.}
\end{minipage} }](img102.png){kind=small}
For this exercise, you will use the last step of the equilibration of ubiquitin in a water sphere, file ubq_ws_eq.restart.coor. You will use the temperature coupling feature of NAMD to set the temperature of the molecules in the outer layer of the sphere to 200 K, while the rest of the bubble is set to a temperature of 300 K. In this way, you will determine the thermal diffusivity by monitoring change of the system's temperature and comparing it to the theoretical expression.
![\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
...ft[-\left(\frac{n\pi}{R}\right)^2Dt\right].
\end{displaymath}}
\end{minipage} }](img103.png)
| cd ../2-6-heat-diff/
|
In order to use the temperature coupling feature of NAMD, you need to create a file which marks the atoms subject to the temperature coupling. The following lines
| tCouple on | |
| tCoupleTemp 200 | |
| tCoupleFile ubq_shell.pdb | |
| tCoupleCol B
|
| cd ../2-6-heat-diff/
|
| mol new ../common/ubq_ws.psf | |
| mol addfile ../1-2-sphere/ubq_ws_eq.restart.coor
|
| set selALL [atomselect top all]
|
| set center [measure center $selALL weight mass]
|
, 
and 
coordinates of the system's center and place their values in the variables xmass, ymass, and zmass:
| foreach {xmass ymass zmass} $center { break }
|
| set shellSel [atomselect top "not ( sqr(x-$xmass) | |
| + sqr(y-$ymass) + sqr(z-$zmass) <= sqr(22) ) "]
|
| $shellSel set beta 1.00
|
| $selALL writepdb ubq_shell.pdb
|
| copy example-output\ubq_cooling.log .
|
| source ../2-3-energies/namdstats.tcl | |
| data_time TEMP ubq_cooling.log
|
Open
, navigate to the directory 2-6-heat-diff/, and choose the file TEMP.dat. Make sure the dropdown menu Files of type is set to All Files.
Chart, choosing XY (Scatter) as your Chart Type, selecting the last Chart sub-type in the lower-right, and clicking Next.
| = 200 | |
| +66.87*(EXP(-0.0146*E$1*C1) +0.25*EXP(-0.25*0.0146*E$1*C1) | |
| +1/9*EXP(-1/9*0.0146*E$1*C1) +1/16*EXP(-1/16*0.0146*E$1*C1) | |
| +1/25*EXP(-1/25*0.0146*E$1*C1) +1/36*EXP(-1/36*0.0146*E$1*C1) | |
| +1/49*EXP(-1/49*0.0146*E$1*C1) +1/64*EXP(-1/64*0.0146*E$1*C1) | |
| +1/81*EXP(-1/81*0.0146*E$1*C1) +1/100*EXP(-1/100*0.0146*E$1*C1))
|
. Click the Series tab and add a series with Y Values: of Column D, Rows 1-201. and X Values: of Column C. Click OK.
You may want to change the line weight of your fit by right-clicking the line in your histogram and selecting Format Data Series.
We have our theoretical temperature expression plotted, but the fit is not optimal. If we change thermal diffusivity in cell E1, we will see the fit change. We will now determine the optimal parameters.
2. Then click and hold the lower-right corner of cell F1 and drag down to row 201 to fill in the error values for each bin.
Solver. If that option is not present, you will need to add in the Solver tool by clicking Tools Add Ins and selecting Solver Add-in from the list.
| Set Target Cell: $G$1 | |
| Equal To: Min | |
| By Changing Cells: $E$1
|
![\begin{figure}\begin{center}
\par
\par
\latex{
\includegraphics[scale=0.5]{pictures/tut_unit02_coolwin}
}
\end{center}
\end{figure}](img105.png) |
s
. You should get a value of around
0.45
cms. How does your result compare
to the thermal diffusivity of water
cms?
, by the following
formula:
, where 
is the specific heat (that was
calculated in Section 2.1.5) and 
is the thermal diffusivity
(assume that the density of the system, 
, is 1 g/mL).
![\framebox[\textwidth]{
\begin{minipage}[r]{0.9\textwidth}
\noindent{\textbf{Ob...
...teins by generating
temperature echoes using MD simulations.}
\end{minipage} }](img113.png){kind=small}
The motions of atoms in globular proteins (e.g., ubiquitin), referred to as internal dynamics, comprise a wide range of time scales, from high frequency vibrations about their equilibrium positions with periods of several femtoseconds to slow collective motions which require seconds or more, leading to deformations of the entire protein.
The internal dynamics of these proteins on a picosecond time scale
(high frequency) can be described as a collection of weakly
interacting harmonic oscillators referred to as normal modes. Since
normal modes are formed by linear superposition of a large number of
individual atomic oscillations, it is not surprising that the internal
dynamics of proteins on this time scale has a delocalized character
throughout the protein. The situation is similar to the lattice
vibrations (phonons) in a crystalline solid. Experimentally, there
exist ways to synchronize, through a suitable signal or perturbation,
these normal modes, forcing the system in a so-called (phase) coherent
state, in which normal modes oscillate in phase. The degree of
coherence of the system can be probed with a second signal which
through interference with the coherent normal modes may lead to
resonances, referred to as echoes, which can be detected
experimentally. However, the coherence of atomic motions in proteins
decays through non-linear contributions to forces between atoms. This
decay develops on a time scale 
which can be
probed, e.g., by means of temperature echoes, and can be described by
employing MD simulations.
In a temperature echo the coherence of the system is probed by
reassigning the same atomic velocities the system had at an earlier
time and then looking to an echo in the temperature at time 
,
as a result of such reassignment. An example is shown in
Fig. 14. At time 
the velocities of all
atoms in the system are quenched; then, at 
the atomic
velocities are reassigned again (quenched) to their values at time
. As a result, a temperature echo, i.e., a sharp dip in
, is detected at a subsequent time 
after 
(i.e., at time 
).
In this section, you will employ MD simulations to generate temperature echoes in ubiquitin by applying the velocity reassignment just described. By modeling ubiquitin as a large collection of weakly interacting harmonic oscillators (normal modes), you will find that:
of the oscillators is about one picosecond;
of the echoes decay exponentially
with the delay time , i.e.,
, the exponential decay being
determined by the dephasing time .
First you will generate temperature echoes in MD simulations, and then you will analyze the echoes in the framework of the normal mode approximation.
The phenomenon of temperature echoes has a simple explanation
within the normal modes description of a protein.
The first velocity reassignment enforces phase coherence for the
oscillators. Because of the frequency dispersion of the normal modes
(as well as the deviation from the harmonic approximation) the phase
coherence of the protein will decay in time
with a characteristic dephasing time
ps.
The second velocity reassignment after a delay time 
is a
``probing signal'' which will test the degree of coherence of the
system at the instant of time it was applied. The depth of the echo
and the instant of time at which it occurs are quantitative
characteristics of the coherence of the internal dynamics of
proteins.
In order to generate temperature echoes, one needs to equilibrate
ubiquitin at the desired initial temperature 
K, e.g., by using
the velocity rescaling method. Once the system is equilibrated, the
thermostat is removed and all the following simulations are carried out
in the microcanonical (NVE) ensemble.
In this section, you will consider two types of temperature echoes. In the first, you will reassign all atomic velocities to zero, and in the second, you will reassign the atomic velocities to a 300 K distribution (the same as the initial temperature of the system) and the second reassignment will be exactly the same as the first.
The temperature quench echo is obtained by resetting to zero
the velocity of all atoms of the protein at both times and
. You should start from the pre-equilibrated protein in
vacuum at 
; the required files are located in the common directory.
You need to run the simulations in the microcanonical ensemble (NVE) by using the configuration file equil.conf in the 2-7-echoes/01_equil_NVE directory. The simulation will run for 500 time steps (fs).
Let's take a look at some features in the equil.conf configuration file. You can read the file using WordPad.
.
| # Continuing a job from the restart files | |
| set inputname ../common/ubi_equil | |
| binCoordinates $inputname.restart.coor | |
| extendedSystem $inputname.xsc
|
| #temperature $temperature | |
| binvelocities $inputname.restart.vel
|
| copy example-output\equil.log .
|
Programs VMD.
Tk Console), navigate to the 2-7-echoes/01_equil_NVE directory if you are not already there.
| source ../../2-3-energies/namdstats.tcl | |
| data_time TEMP equil.log
|
| source auto-corr.tcl
|
Open, navigate to the directory 2-7-echoes/01_equil_NVE/, and choose the file auto-corr.dat. Make sure the dropdown menu Files of type is set to All Files.
Chart, choosing XY (Scatter) as your Chart Type, selecting the last Chart sub-type in the lower-right, and clicking Next. A preview should show that Excel assigned Column A to x-axis and Column B to the y-axis, as we would like. Click Finish.
You will now fit the simulated temperature autocorrelation function to a decaying exponential approximation (see the Science Box).
| =EXP(-A1/D$1)
|
, and we will alter it to obtain a better fit. Enter an initial guess of 1 in cell D1. Then extend the formula in C1 down to cell C26.
. Click the Series tab and add a new series with Y Values: of Column C, Rows 1-26. and X Values: of Column A. Click OK.
You may want to change the line weight of your functions by right-clicking the lines in your graph and selecting Format Data Series.
We have our approximation to the temperature autocorrelation function plotted, but the fit is not optimal. If we change decay time in cell D1, we will see the fit change. We will now determine the optimal decay time.
2. Then click and hold the lower-right corner of cell E1 and drag down to row 26 to fill in the error values for each bin.
Solver. If that option is not present, you will need to add in the Solver tool by clicking Tools Add Ins and selecting Solver Add-in from the list.
| Set Target Cell: $F$1 | |
| Equal To: Min | |
| By Changing Cells: $D$1
|
![\begin{figure}\begin{center}
\par
\par
\latex{
\includegraphics[scale=0.5]{pictures/tut_autotempwin}
}
\end{center}
\end{figure}](img134.png) |
| cd ..\02_quencha
|
. You are going to quench the temperature (set it to zero), and monitor the recovery of the system as it runs for time steps. If you would like to probe other values of , you should edit your configuration file echo.conf using WordPad, changing the value of tau:
| set tau 200
|
| copy example-output\echo.log .
|
| data_time TEMP echo.log
|
Now, you will quench the system for the second time.
Note that the configuration file echo.conf indicates that the simulation will run for 
steps: run [expr 3*$tau]
| copy example-output\echo.log .
|
| data_time TEMP echo.log
|
| move TEMP.dat temp3.dat
|
| copy ..\01_equil_NVE\TEMP.dat .\temp1.dat | |
| copy ..\02_quencha\TEMP.dat .\temp2.dat
|
| type temp1.dat > temp.dat | |
| type temp2.dat » temp.dat | |
| type temp3.dat » temp.dat
|
You now have the temperature data for the whole temperature quench experiment. You will use Excel to analyze the data.
Open, navigate to the directory 2-7-echoes/03_quenchb/, and choose the file temp.dat. Make sure the dropdown menu Files of type is set to All Files.
Chart, choosing XY (Scatter) as your Chart Type, selecting the last Chart sub-type in the lower-right, and clicking Next. A preview should show that Excel assigned Column A to x-axis and Column B to the y-axis, as we would like. Click Finish.
Note anything unusual? You quenched the temperature twice, and yet you see three drops in temperature! The third one is a temperature echo. See Fig. 16.
Next, you will compare in the vicinity of the echo () obtained from the MD simulation with the theoretical prediction (9) involving the temperature autocorrelation function obtained also from the MD simulations. The degree of agreement between these two results is a measure of the accuracy of the harmonic approximation.
| source harmonic.tcl
|
Get External Data Import Text File, navigate to the directory 2-7-echoes/03_quenchb/, and choose the file harmonic.dat. Make sure the dropdown menu Files of type is set to All Files.
. Click the Series tab and add a series with X Values: of Column C, Rows 1-601. and Y Values: of Column D, Rows 1-601. Click OK.
| = 75 - MIN(B800:B1000)
|
=300K. Note that this formula is designed knowing the echo is between 800 and 1000 fs.
![\begin{figure}\begin{center}
\par
\par
\latex{
\includegraphics[scale=0.5]{pictures/tut_quenchechowin}
}
\end{center}
\end{figure}](img140.png) |
You have calculated the echo depth for a particular value of the delay time . You may wish to repeat the experiment using other values of (100-800 fs). (You must alter the harmonic.tcl script.) If this is done and the data is plotted and fitted to a single exponential
, one obtains the so-called dephasing time 
, which is an inherent property of the system.
Here you will generate a temperature echo by replacing the
velocities of the atoms at time with the ones you assign to
them at time 
. Moreover, by choosing these velocities according
to the Maxwell-Boltzmann distribution corresponding to K,
i.e., the temperature of the equilibrated system (
), there will
be no discontinuities (jumps) in at and , yet you
will observe a temperature echo at a later time, i.e., a sharp dip in
the temperature at time 
. Your goal is to determine
and the depth 
of the echo. To this end, you will
need to follow a similar procedure to the one in the previous
exercise.
according to a
Maxwell-Boltzmann distribution for a temperature 
, saving
the reassigned velocities to a file 300.vel.
In order to do this, the following is added to the end of your
configuration file echo.conf:
| run 0 | |
| output 300
|
| run $tau
|
| copy example-output\echo.log .
|
| 500 300.5656 | |
| 500 300.5656 | |
| 501 301.0253 | |
| 502 302.5395 | |
| ...
|
| velocities ../04_consta/300.vel
|
time steps.
| copy example-output\echo.log .
|
| move TEMP.dat temp2.dat | |
| copy ..\04_consta\TEMP.dat .\temp1.dat
|
| type temp1.dat > temp.dat | |
| type temp2.dat » temp.dat
|
| Echo depth = 300 - Minimum Temperature
|
fs. How would you explain your findings?
![\begin{figure}\begin{center}
\par
\par
\latex{
\includegraphics[scale=0.5]{pictures/tut_quench-schem}
}
\end{center}
\end{figure}](img124.png)
![\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
...0)]>^2}
\nonumber \\ &\approx &exp(-t/\tau_0)
\end{eqnarray*}}
\end{minipage} }](img132.png)
![% latex2html id marker 13272
\fbox{
\begin{minipage}{.2\textwidth}
\includegr...
...{1}{2}C_{T,T}(\vert t-2\tau\vert)\right]\;.
\end{equation} }
\end{minipage} }](img138.png)