0README

atom-fluctuations
|-- 0NOTE                obvious
|-- 0README              obvious
|-- yout184-dget.Data    trajectory data:  chi-2 of BPTI Tyr-35
|-- yout184.s            R commands for analysis of trajectory
|-- yout184-out.s        same, but to generate ps and txt files
|-- yout184-out.ps       plots generated by analsysis
`-- yout184-out.txt      text generated by analsysis



TRAJECTORY AND MODEL FOR ANALYSIS

A molecular dynamics simulation for a protein at equilibrium generates
a trajectory for a set of variables sufficient to describe the motions
of the molecule.  These stochastic motions, driven by the ambient
thermal energy, are fluctuations about the equilibrium values.  The
set of variables can be the atom coordinates, from which at each
time point one can calculate trajectories for particular bond
lengths, bond angles, and dihedral angles.

We consider fluctuation of the dihedral angle chi-2 of tyrosine 35
of bovine pancreatic trypsin inhibitor (BPTI).  These motions
correpond to rotations of the tyrosine ring about an axis approximately
coincident with the C-beta--C-gamma bond.  Tyr-35 is largely buried
in the protein matrix.  We can view it as moving in a potential
well the shape of which is determined by interactions of the ring
atoms with atoms of residues adjacent in tbe protein structure.
The thermal motions of the adjacent atoms provide both the random
force driving the ring motions and the viscous force opposing its
motion (relaxation).

This picture suggests that the fluctuations of chi-2 of Tyr-35 may
be modeled as those of an harmonically bound Brownian particle.
The harmonic motion, described by a torsional force constant, or
equivalently, a frequency of oscillation, is determined by the shape
of the potential well.  The Brownian motion of the ring is determined
by thermal motions of the protein matrix and is described by a
friction constant, in effect modeling the protein matrix as a solvent
in which the ring is embedded.

The trajectory of chi-2 of Tyr-35 of BPTI analyzed below is from a
simulation at 100K (please do not ask, why not 300K?).  Details of
the simulation are in file "0NOTE".  Part of the analysis is fit
of the data with a tethered Brownian oscillator model.  The quality
of the fit is adequate.  It is perhaps remarkable that a simple
2-parameter model can capture the essential features of Brownian
motion of a ring buried in the protein matrix.


R LANGUAGE COMPUTATIONS AND RESULTS

The file "yout184.s" has a set of commands to analyze a trajectory
for the dihedral angle chi-2 of Tyr-35 of BPTI:

	read in the chi-2 trajectory as a vector of angular
	fluctuations (degrees from the mean value) for each time
	point, saved in an S/R dump format (dput/dget);

	construct an equal-length vector of the time of simulation
	for each time point (picoseconds from start of run);

	two plots of the trajectory:

		full trajectory, 16384 data points, 32.768 ps;

		first 500 points, 1 ps;

	Fourier transform the signal and normalize;
	two plots, of modulus = Mod(fft(yout184)/16484)
		
		full transform, 8192 frequency bins, F-hat and F-hat*;
			(2nd half of plot is redundant)

		first 750 points of the transform, F-hat;

	construct the time correlation function:

   C(t) = Re(fft(Conj(fft(yout184)/16484)*fft(yout184)/16484,inverse=TRUE))

	fit the tethered Brownian oscillator model to the data for
	C(t) (first 500 points), using the nonlinear least squares
	algorithm of nls(); the data are from a simulation for 100 K:


	mean chi-2 =      -3.845215e-05
	variance chi-2 =  19.49047
	std.dev. chi-2 =  4.414802

	Parvseval's theorem - size signal vector invariant to change basis
	var(yout184) = sum(yout184^2)/16384 = sum(Mod(a)^2)) = C(0) =  19.48929

	first 20 values of C(t); note C(0) = variance
 [1] 19.489285 18.009616 15.247437 13.633616 13.706762 14.740587 16.006862
 [8] 16.270025 14.439963 11.713520 10.733692 12.100012 13.886611 14.502295
[15] 13.789323 11.949852  9.803482  9.191256 10.793810 12.592394


	summary of fit to C(t)

	Formula: Ct ~ C * exp(-time/tau) * (cos(sqrt(XX/C - 1/tau^2) * time) + 
	    1/(tau * sqrt(XX/C - 1/tau^2)) * sin(sqrt(XX/C - 1/tau^2) * 
		time))

	Parameters:
	     Estimate Std. Error t value Pr(>|t|)    
	tau  0.185765   0.006422   28.92   <2e-16 ***
	C   12.316130   0.184796   66.65   <2e-16 ***
	---
	Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

	Residual standard error: 1.664 on 498 degrees of freedom

	Correlation of Parameter Estimates:
	     tau
	C 0.2501


	torsional spring constant, from fit value C(0): K  =  52.77597  Kcal/mol

	freq. of torsional oscillation, from K: w0 =  16.93435 ps^-1

	freq. of Brownian model oscillation, from C & tau: w1 =  16.05596 ps^-1


ANALYSIS

Plots 1 and 2 of file "yout184-out.ps" and the averages listed first
in file "yout184-out.txt" (see above values) describe the trajectory
of chi-2 of BPTI, from a simulation at 100K.  The tyrosine ring
undergoes Brownian motion with an RMS deviation of 4.4 degrees.
Inspection of plot 2 indicates the relaxation time for a large
fluctuation is between .1 and .2 ps.

The Fourier transform of the signal is shown in plots 3 and 4.  The
strongest contributions are at the frequencies: 130, 75, 30, 16,
and 7 ps^-1 (corresponding to 4300, 2500, 1000, 525, and 230 cm^-1).
Below 4 ps^-1, the transform shows white noise (all frequencies
with the same amplitude), consistent with Brownian motion of the
ring, driven by very fast random collisions of the ring with its
environment (the time scale of the collisional events is much shorter
and well separated from that of the slow Brownian diffusion, so
that ring motions can be modeled as driven by a delta-function
random force).

The time correlation function C(t) is obtained directly from the Fourier
transform, as the inverse transform of the convolute:

	C(t) = F^-1(F(chi-2(t))* * F(chi-2(t)))

where F() is the Fourier transform operator, F()* its complex
conjugate, and F^-1() the inverse.

The tethered Brownian oscillator model was fit to C(t):

  C(t) = C(0) * exp( -t/tau) * ( cos( w1 * t ) + 1/( tau * w1 ) * sin( w1 * t )

where w1 = sqrt(K/I - 1/tau^2), K=kT/C(0), and I is the moment of
inertia of the ring.  The two free parameters are tau, a relaxation
time, and C(0), the zero-time correlation function.

Plot 5 shows C(t) vs time, with the dashed line calculated from the
best-fit parameters (tau=0.186, C(0)=12.3).

The fit is good up to t ~ 0.6 ps; for longer times the motion is
less strongly damped than predicted by the model.  The fit value
of C(0) is one third less than the value calculated from the
trajectory.  We may view the model as describing the principal
features of the ring motion, while recognizing that the motion has
additional components.  As noted, it is remarkable that a simple
two-parameter model can account for the essential features of so
complex a system as a tyrosine ring moving within an asymmetric and
dynamic protein matrix.

The relaxation time tau is consistent with the decay of large
fluctuations seen in the chi-2 trajectory (plot 2).

From the estimate for the parameter C(0) of the fit (12.3 degrees^2),
we obtain the torsional force constant for oscillation of the
tyrosine ring:

	K = RT/1000*1/C(0) = 53 Kcal/mol

which corresponds to a frequency of oscillation (calculated for a
moment of inertia I = 7.7 10^15 gm-cm^2),

	w0 = sqrt(K/I) = 16.9 ps^-1

This oscillator frequency is close to the Fourier component of 16 ps^-1.

The frequency, w1 = 16.1 ps^-1, of the periodic term of the tethered oscillator
model is approximately equal to w0, i.e., 1/tau^2 << w0.

The torsional force constant, K,  corresponds to a soft motion, of
much lower frequency than for bond stretching or bond bending,
reflecting the deformability of the protein matrix and the weakness
of non-covalent interactions.


John Rupley
 rupley@u.arizona.edu -or- jar@rupley.com
 30 Calle Belleza, Tucson AZ 85716 - (520) 325-4533; fax - (520) 325-4991
 Dept. Biochemistry & Molecular Biophysics, Univ. Arizona, Tucson AZ 85721
