| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
Biophys J, October 2002, p. 1854-1866, Vol. 83, No. 4
Biochemistry, Molecular Biology, and Biophysics Department, University of Minnesota, Minneapolis, Minnesota 55455 USA
 |
ABSTRACT |
|---|
 TOP
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
RESULTS AND DISCUSSION
CONCLUSION
REFERENCES
|
|---|
We have developed a computational molecular dynamics technique to simulate the motions of spin labels bound to the regulatory domain of scallop myosin. These calculations were then directly compared with site-directed spin labeling experimental results obtained by preparing seven single-cysteine mutants of the smooth muscle (chicken gizzard) myosin regulatory light chain and performing electron paramagnetic resonance experiments on these spin-labeled regulatory light chains in functional scallop muscle fibers. We determined molecular dynamics simulation conditions necessary for obtaining a convergent orientational trajectory of the spin label, and from these trajectories we then calculated correlation times, orientational distributions, and order parameters. Simulated order parameters closely match those determined experimentally, validating our molecular dynamics modeling technique, and demonstrating our ability to predict preferred sites for labeling by computer simulation. In several cases, more than one rotational mode was observed within the 14-ns trajectory, suggesting that the spin label samples several local energy minima. This study uses molecular dynamics simulations of an experimental system to explore and enhance the site-directed spin labeling technique.
| |
INTRODUCTION |
|---|
|
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS AND DISCUSSION
CONCLUSION
REFERENCES
|
|---|
Site-directed spin labeling (SDSL) is a powerful
technique to probe protein dynamics and conformation (Hubbell and
Altenbach, 1994
; Hubbell et al., 2000). Using site-directed cysteine
mutagenesis, spin labels are attached throughout a protein, one site
per sample, and electron paramagnetic resonance (EPR) spectra are used
to measure spin label rotational mobility and accessibility and to interpret this data in terms of protein secondary and tertiary structure. There is a large body of evidence to show that the restriction of nanosecond dynamic disorder exhibited by a spin label
attached to a protein is a direct reflection of the local structure and
dynamics of the protein, due to nearby residues within the protein or
to contacts made with other molecules. Indeed, the size of a binding
pocket can be measured by probe mobility as the steric effect of
differently sized substituents (Columbus et al., 2001).
When attached to a specific protein residue, probe restriction should
be well modeled by a molecular dynamics (MD) simulation reflecting the
restrictions imposed by neighboring protein atoms. EPR spectroscopy has
the right spectral resolution and time scale to test the validity of MD
simulations. The extent of narrowing of the nitroxide EPR spectrum is
optimally sensitive to the amplitude (measured by the order parameter)
of nanosecond and subnanosecond rotational motions. It should be
possible to calculate the order parameter accurately from an MD
simulation, provided that the probe orientational trajectory reaches
equilibrium within a few nanoseconds. Steinhoff and Hubbell have shown
that MD trajectories of a nitroxide spin label attached to a tripeptide
with fixed C
s generally reach equilibrium
within a few nanoseconds and yield EPR simulations that show good
qualitative agreement with whole-protein EPR experiments in which
structurally homologous residues were labeled (Steinhoff and Hubbell,
1996). In the present study we apply a similar approach to MD
simulations based directly on a protein crystal structure.
Spin-labeled muscle fibers as a model system
Muscle fibers provide an excellent system in which to test both
the SDSL technique and the accuracy of modeling using MD. When spin
labels are rigidly and stereospecifically bound to sites within the
muscle fiber, experimentalists can take advantage of the highly
ordered, oriented system to study rigid-body motion of protein
structures (Thomas and Cooke, 1980). In the case of an unoriented
sample, measurements of probe mobility can directly reflect local
conformational changes within muscle proteins (Ostap et al., 1993). The
regulatory light chain (RLC) of myosin can be easily exchanged in a
muscle fiber, offering the opportunity to use molecular biology to
engineer specific spin-labeling sites to investigate local structure or
global movement (Baker et al., 1998; Roopnarine et al., 1998).
EPR spectroscopy of spin-labeled myosin has demonstrated that the
regulatory domain undergoes a ~36° rotation upon muscle contraction, acting as a lever arm to produce mechanical force. This
result was found by spin labeling the single native cysteine residue
(Cys-108) of the chicken gizzard RLC incorporated into myosin in
scallop muscle fibers (Baker et al., 1998). The spin label
3-(5-fluoro-2,4-dinitroanilino)-2,2,5,5-tetramethyl-1-pyrrolidinyloxy (FDNASL) binds specifically and rigidly to Cys-108, resulting in a
probe that is sensitive to the orientation of the regulatory domain
relative to the muscle fiber. We sought to complement this SDSL study
by systematically making several site-directed single-cysteine mutants
on the RLC, attaching FDNASL, and measuring their EPR spectra in an
oriented muscle fiber, thus potentially mapping out orientational
conformations of the RLC. However, it remained a significant challenge
to choose labeling sites that would yield restricted probes. There is a
low probability that any one Cys mutant will possess both functionality
and a rigidly bound, oriented probe, so identifying suitable residues
by trial-and-error mutagenesis, expression, and purification of
individual RLC mutants is very time consuming.
We reasoned that a MD simulation of a spin-labeled protein might
shorten this process by predicting which labeling sites will produce
the desired level of orientation and mobility. By performing these
simulations on an experimentally accessible system, the present study
offers a direct test of the MD calculation, representing a significant
advance over previous studies of spin label dynamics (Steinhoff and
Hubbell, 1996). Order parameters of specific spin-labeled mutants can
be predicted and then compared directly with order parameters measured
experimentally, verifying both the SDSL and MD techniques.
A detailed MD simulation of spin labels also presents us with ways to
supplement EPR spectral measurements to resolve otherwise ambiguous
conclusions. Absolute orientation of probe trajectories with respect to
a known reference frame can give information about the conformations of
the proteins to which they are attached. It is also possible to
simulate EPR spectra directly from probe trajectories and compare these
with measured spectra (Freed, 1976; Robinson et al., 1992; Steinhoff
and Hubbell, 1996). This presents the possibility of resolving
ambiguous spectral shapes with unprecedented detail, e.g., by
separating the competing effects of orientational disorder and
fast-motional EPR tensor averaging. Furthermore, by performing long MD
simulations at many different labeling sites, we can begin to get a
picture of the general properties of attached spin probes and what
factors are important in determining their dynamic behavior.
In the present paper, we first describe our computational methods and then examine in detail the dynamic behavior of a series of spin-labeled myosin RLC mutants over long MD simulations (14 ns). MD simulation is shown to accurately predict observed order parameters as measured by EPR spectroscopy, providing an excellent way to predict the usefulness of proposed site-directed mutants.
| |
MATERIALS AND METHODS |
|---|
|
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS AND DISCUSSION
CONCLUSION
REFERENCES
|
|---|
Experimental system
Sample preparation
RLC was purified from fresh chicken gizzard myosin (Persechini and Hartshorne, 1983), then spin-labeled with FDNASL as described previously (Baker et al., 1998). Muscle fiber bundles were prepared from Placopecten magellanicus (sea scallop) adductor muscle,
and the native RLC was exchanged for the labeled chicken gizzard RLC, as described previously (Brust-Mascher et al., 1999). These fiber bundle samples were placed in capillaries for EPR spectral acquisition. Minced fibers were prepared by cutting labeled fibers into fragments with dissection scissors and were then flowed into capillaries and
trapped with glass wool or placed in the well of a flat cell for
spectral acquisition.
RLC mutant expression and purification
Chicken gizzard RLC mutants were prepared using a null cysteine (C108A) cDNA construct, kindly provided by Dr. Renne Lu from the Boston Biomedical Research Institute. This cDNA was ligated into the pet3a expression system (Novagen, Madison, WI) between the sites NEI and EcoRI. Using the QuikChange Mutagenesis system (Stratagene, La Jolla, CA), single cysteines were mutated at desired positions in the RLC. The pet3a-RLC with the desired mutation was then cotransformed with the pLys plasmid (Novagen) into the Escherichia coli expression strain, JM109DE3. The pLys plasmid encodes for T7 lysozyme, which shuts down basal transcription and stabilizes the pet3a-RLC plasmid. For each mutant, a 10-mL overnight culture was used to inoculate a large 1-L culture. The cells were incubated, shaking at 37°C until an optical density of 0.6 was reached, and then induced with IPTG (isopropyl 1-thio-
-D-galactopyranoside). The culture was
incubated at 37°C for another 3 h, and the cells were harvested
by centrifugation at 5000 × g for 10 min. The cells were resuspended in 50 mL of 25 mM Tris-HCl, pH 8.0, 5 mM EDTA, 50 mM
glucose, 1 mM phenylmethylsulfonyl fluoride, and 0.2 mg of lysozyme.
The suspension was incubated on ice for 1 h and subjected to one
cycle of freeze/thaw using an isopropanol/dry ice bath. The cells were
then incubated with 10 mM MgCl2 and 5 µg/mL
DNase for 1 h. Triton was added to 0.1%, and the cells were
centrifuged at 15,000 × g for 15 min. Two more washes
with the above buffer and 0.1% Triton were carried out, and one final
wash was done without Triton. The pellet was dissolved in 10 mL of 7 M
urea, 30 mM Tris-HCl, pH 7.5, 50 mM NaCl, and 1 mM dithiothreitol and centrifuged at 100,000 × g for 30 min. A DE-52
(Whatman, Clifton, NJ) anion exchange column was equilibrated with 4 M
urea, 30 mM Tri-HCl, pH 7.5, 50 mM NaCl, and 1 mM dithiothreitol. The
supernatant was applied to the column, and a linear gradient from 0.05 to 0.3 M NaCl was used to elute the RLC. The fractions containing RLC
were pooled and dialyzed against 5 mM ammonium bicarbonate and then lyophilized.
EPR spectroscopy
X-band EPR spectra were acquired on a Bruker ESP 300 spectrometer or a Bruker EleXsys 500 spectrometer (Bruker Instruments, Billerica, MA). Spectra of oriented muscle fibers and randomly oriented minced fibers were acquired as described previously. Spectra were obtained in rigor buffer (20 mM MOPS, 5 mM MgCl2, 1 mM EGTA, and 0.1 mM NaN3, pH 7).Molecular dynamics model
Base structure and molecular model
Because a high-resolution crystal structure of chicken gizzard regulatory domain is unavailable, the known sequence (Maita et al., 1984) was homology-modeled onto the 2.8-Å resolution crystal structure
of scallop myosin regulatory domain (Xie et al., 1994, PDB structure
1scm) using the homology module of InsightII/Biosym. Sequence homology between the scallop regulatory light chain (Maita et
al., 1984) and chicken gizzard regulatory light chain (Maita et al.,
1981) was determined by BLAST to be 47% identical and 70% similar.
The final molecular model consists of 5597 atoms: 359 residues
including a section of the heavy chain (M773-L837), the ELC (L3-P154),
the RLC (L12-E153), and two Ca2+ ions. Disordered
residues not in the crystal structure were not included in the model.
Attaching the probe to the structure
To model a particular spin-labeled mutant, the native residue was replaced by cysteine, the native cysteine replaced by alanine, and the probe molecule covalently attached to its associated sulfur atom. In each case, the mutated residue was appended to the C terminus of the RLC, which conserves molecular mass across all the mutants so their energies can be compared. Models of the wild-type structure (C108) and several mutants produced by site-directed mutagenesis were constructed: A105C, M129C, and D131C. Other "predicted" mutants, G65C, M72C, and G95C were first constructed computationally and later expressed experimentally (Fig. 1 A).
|
Energy minimization and MD performed locally around the probe
Energy minimization and MD were calculated using the Discover3 engine from InsightII/Biosym. We restricted our attention to local minimization and dynamics around the probe; this is justified because experimental mutants are completely functional and we expect the regulatory domain to be minimally perturbed by the presence of the probe. Due to the truncated nature of the model, both from the missing residues in the crystal structure and its detachment from the larger myosin S1 complex, it is prudent to consider only the local dynamics of the probe. As documented in Results and Discussion, we determined that a simulation including atoms within a 20-Å radius around the probe gave accurate results, after examining convergence criteria from simulations of radii ranging from 14 to 24 Å. Force field parameters. The CVFF force field (Dauber-Osguthorpe et al., 1988) was used for energy minimization and
dynamics. Standard atomic potentials for peptides were assigned to the
homology-modeled structure after correcting new double bonds. Unique
partial charges for the probe were calculated with the program MOPAC,
using Mulliken partial charges for FDNASL (Fig. 1 B) using
the AM1 Hamiltonian (Dewar et al., 1985). To account for the unusual
electronic structure in the nitroxide group, the nitrogen atom was
assigned parameters for a partial double bond, preserving the planar
orientation of the oxygen. All simulations were calculated in vacuo
with no counter ions and no ionized residues.
Energy minimization. For each structure, atoms within 20 Å of the probe were energy-minimized in vacuo using the conjugate gradient method to within 0.01 kcal·mol
1·Å1. The
same method was then applied to the entire molecule for 50,000 iterations, yielding a maximum derivative of ~1.5
kcal·mol1·Å1 for
the entire protein. This accounted for larger-scale relaxations upon
perturbation by a probe.
Molecular dynamics. MD calculations were performed within a
subset of allowed motion, defined to be all atoms within 20 Å from the
nitroxide N in FDNASL. If a Ca2+ ion was within
this distance, its position was fixed throughout the simulation. Atoms
more than 45 Å from the probe were ignored, whereas atoms between 20 and 45 Å were included in the simulation but fixed in position.
Trajectory points were calculated in 3-fs time steps using the Velocity
Verlet method. The RATTLE algorithm was used to constrain hydrogen bond
lengths to a tolerance of 105. Using this, the
simulation was stable with a 3-fs time step. The trajectory of the
entire simulation was written to file every 20 ps. Nonbonded force
field terms were calculated using the cell multipole method with a
10-Å cutoff, 1.0-Å update width, and fixed intrinsic dielectric
constant 
= 1.0. The thermodynamic ensemble was
(constant number, volume, and temperature), T = 298 K ± 10 K tolerance by velocity scaling. All simulations were
run on a network of queued SGI Origin 200 computers, on multiple 195 MHz R10000 processors.
In most cases, 14-ns trajectories were assembled from shorter 1, 2, and
4 ns trajectories, each continued from a "restart" file. At these
"restart" points, the simulation was resumed with the previous atom
positions from the last time step, but with atom velocities reassigned
randomly from a Boltzmann distribution at 300 K. This provided a
convenient way to test the stability of the simulation with respect to perturbation.
The nitroxide EPR spectrum is determined almost entirely by the
time-dependence of the angles 
and 
that define the orientation of the applied magnetic field H with respect to the principal axis
system of the nitroxide group, as indicated in Fig.
2. Therefore, to enable calculation of
EPR parameters from the trajectory, coordinates of the spin probe N-O
vector and the probe's p-orbital z-axis, defined as the
cross-product of the C2-N and
C5-N vectors (Fig. 2), were written to file at
every time step.
|
Alignment of light chain domain
To determine the orientation of the spin label relative to the muscle fiber (actin filament) axis (which is aligned with the magnetic field H in our experiments), we started from the structure of an actin filament decorated with chicken skeletal S1 (Mendelson and Morris, 1997), and aligned our MD model of the spin-labeled regulatory domain
onto the structure using the 19 residues in the heavy chain closest to
the catalytic domain, 782 to 800. Although the length of the heavy
chain in our MD model extends beyond residue 800, we did not use those
residues for alignment, because the structures of the regulatory
domains for skeletal and scallop muscle appear to diverge beyond that
point (Houdusse and Cohen, 1996).
Calculation of order parameters from probe trajectories
To compare quantitatively the results of the MD simulation with experimental EPR data, we calculated the spin label order parameter
|
(1) |
cos2
= 
cos2 
() sin d,
directly from the MD trajectory, in which is the angular deviation
from the mode of the orientational distribution of the z
axis of the nitroxide (Fig. 3). This
order parameter determines the EPR spectrum for rapid restricted
rotational diffusion in a randomly oriented sample, as discussed below.
The distribution () is calculated by first determining the
distribution of values (Fig. 2), then defining = 0 at the mode of this distribution. A Windows C++ program called Dynamics
was used to analyze and display trajectory data. The program also
calculates the probe's angular distribution with respect to the muscle
fiber axis, which corresponds to the distribution of the angles and
in Fig. 2 in the case of a muscle fiber oriented parallel to the
applied field.
|
T'
)
to the maximum possible tensor values (TT) (Gaffney,
1976; Squier and Thomas, 1989):
|
(2) |
|
Autocorrelation of probe trajectories
We define the autocorrelation of a probe trajectory of length T to be
|
(3) |
(t) is the angle of the z axis
of the nitroxide with the mode of the trajectory. We take the
autocorrelation over the cosine of the mode angle because of the useful
quantities G (0)=cos2 and
lim
G(
) = cos2, from which order parameter can be
calculated (Eq. 1). The rotational correlation time is defined from the
best fit to G() = G(0)exp(/c) over a specified
time window.
Spin label and side-chain radial distributions, RMS, and diffusion constants
To understand the relationship between the calculated order parameter and the spin label's environment, we analyzed the distances from the spin label and adjacent residues across the MD trajectory. For groups of atoms a and b, INSIGHT calculates the radial distribution functionab(r), the probability that atoms
from groups a and b are separated by distance
r. Then the rms distance
rrms is calculated from
|
(4) |
by
|
(5) |
|
(6) |
).
The values of rrms and D,
between a residue and a spin label segment, were calculated across 698 time points along a 14-ns trajectory (every 20 ps). Atom groups for
nearby residues were defined as all atoms in the side chain including
the C. Distances were calculated
from each residue to the entire spin label and also to each of its two
major segments.
| |
RESULTS AND DISCUSSION |
|---|
|
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS AND DISCUSSION
CONCLUSION
REFERENCES
|
|---|
Convergence and consistency of MD trajectories
We used C108-FDNASL as a test case to demonstrate that our
simulation produces convergent, self-consistent dynamics trajectories that are insensitive to velocity perturbations (Allen and Tildesley, 1987). We also measured the sensitivity of the simulation to starting conformer.
Most trajectories converged within 1 ns with potential energy and rms differences remaining stable over the entire 14-ns trajectory (Fig. 5), giving equivalent results in independent runs. The potential energy (Fig. 5 A) and RMS differences from the starting structure (Fig. 5 B) converged within 500 ps. With the exception of A105C-FDNASL (discussed below), the trajectories were insensitive to velocity perturbations (Fig. 5 C).
|
Determining an appropriate simulation subset of allowed motion
To determine the appropriate size of the movable subset of atoms, we calculated 14-ns trajectories in which the movable subset had radii ranging from 14 to 24 Å (Fig. 6). All dynamics trajectories had nearly identical average conformations. The principal correlation time, determined by fitting the autocorrelation function to a single exponential, was less than 50 ps in each case. This rapid convergence of the orientational distribution justifies the calculation of an order parameter S (Eq. 1) that characterizes the amplitude of subnanosecond probe rotational motion. The order parameter decreased significantly with increasing radius, suggesting that artifactual motional constraints were significant at short radii but converged for radii of 18 Å and greater (see Fig. 8). Therefore, we used a radius of 20 Å in all subsequent simulations. These results highlight the importance of examining the effects of simulation constraints to avoid significant artifacts.
|
Convergence sensitivity to starting conformation
C108-FDNASL was used as a test case to examine the sensitivity of
convergence to the starting conformation. When the initial position of
the spin label was manually repositioned in the minimized structure to
five different conformations, all trajectories converged within 500 ps
to essentially identical structures and energies. To examine
convergence sensitivity with more distant starting conformations, a
200-ps trajectory of C108-FDNASL at 600-K was simulated, and eight
low-potential energy conformers were chosen from approximately every 50 ps in the trajectory. This created a series of conformations that were
increasingly farther away on the energy landscape from the
original minimized starting structure. These conformations were then
used as starting structures for subsequent dynamics for 6 ns at 300 K. The trajectories uniformly converge over the first 1 ns to single-mode
dynamics about several local minima stable over the 6-ns simulation.
The three lowest-energy trajectories that were computed
each with an
RMS difference from the starting structure of less than 1.3 Åreproduced the trajectory originally determined by our minimization
and dynamics methods, whereas the higher-energy trajectories did not.
Besides having higher potential energies, the trajectories about more
distant local minima were distinguished by their lower calculated order parameters across the last 2 ns of the trajectory (Fig.
7) and distortion of the peptide backbone
conformation incongruent with the crystal structure. No lower-energy
minima were found. We conclude that C108-FDNASL has a well-defined
energy minimum determined by our methods and that starting
conformations different by more than an RMS distance of ~1.3 Å from
the average structure will not produce the correct dynamics within the
length of our simulation.
|
Trajectories of FDNASL-RLC mutants
Shown in Fig. 8, for each spin-labeled site, is the experimentally determined EPR spectrum, the simulated 14-ns MD trajectory, and the angular distribution of the probe.
|
Commonalities among the trajectories
The trajectories (Fig. 8, center) and their corresponding
distribution functions (Fig. 8, right) reveal the number of local orientational states and the extent of rotational disorder about each
state. In several cases, the probe orientation shows evidence of
jumping between discrete orientational states, whereas the energy of
the entire system shows a single well-defined energy minimum.
Inspection of the structures reveals a simple explanation: the primary
source of motion among the reorientational states is the rotation of
the relatively rigid five-membered pyrrolidinyl ring about its proximal
amine bond. This has been suggested previously for other spin labels,
based on MD simulations (Steinhoff and Hubbell, 1996) and experiments
(Columbus et al., 2001; Langen et al., 2000). Whereas internal steric
hindrance in the spin label favors a trans conformer of the
pyrrolidinyl ring with respect to the amine group (as shown in Fig. 1
B), thermal motion is apparently enough to overcome this
activation barrier in favorable cases.
We will now discuss briefly the distinguishing characteristics of each spin-labeled site.
C108
EPR experiments have shown that the spin label attached to this site is quite restricted in its rotational mobility with respect to the protein (Baker et al., 1998) with an order parameter of 0.82. The
simulated MD trajectory displays similar restricted probe dynamics,
approximating a single Gaussian mode (Fig. 8). Fourteen nanoseconds of
MD with velocity perturbation failed to produce dynamics about any
other local minima. Inspection of the structure shows that the probe
rests in a pocket between the two light chains and is stabilized by
steric interactions with two nearby leucines on the heavy chain and
possible van der Waals interactions between the phenyl ring of FDNASL
and nearby F109.
D131C
The spin label attached to C131 resides in a relatively exposed pocket between the two light chains. As a result, long bimodal fluctuations of the pyrrolidinyl ring are seen on the order of ~1 ns. Approximately 9 ns into the trajectory, nearby R813 on the heavy chain assumes a new rotameric conformation, causing a small change in steric boundaries, and a slight change in the conformation of the entire spin label in the last 5 ns of the trajectory, while maintaining the bimodal rotations of the nitroxide ring. The two average spin label atom positions are shown in Fig. 9 A. The shift can also be seen in the trajectory (Fig. 8) as a sudden increase of ~0.5 radians in the angle at 9 ns. In the actual protein, these shifts probably occur as side chains exchange between different rotamers, partially blurring the distinct modes. This
is apparent in the experimental EPR spectrum, which may display motional averaging.
|
M60C
Like D131C-FDNASL, the trajectory of the spin label attached to C60 displays bimodal dynamics but on a much faster timescale on the order of ~81 ps (Fig. 8; Table 1). The spin label's partially exposed conformation on the RLC and fast bimodal time scale indicate that its nearby groups assume a conformation that lowers the energy barrier between rotameric states of the pyrrolidinyl ring.
|
M72C
The spin label FDNASL bound to C72 shows fluctuating orientation and slow energy equilibration over the first 4 ns of this trajectory. During this equilibration period, the spin label appears to explore two orientations inside a crowded binding pocket. The spin label displays dynamics about only one of these orientations for the rest of the simulation, but the relatively high average potential energy suggests that there is a more stable global energy minimum that has not been reached.M129C
This probe lies between the essential and regulatory light chains and is partially exposed. The MD displays single-mode dynamics but with a larger amplitude than at C108. A tryptophan residue (W21) in the essential light chain is the closest to the nitroxide group and has an average conformation showing a roughly planar configuration with the methyl groups on FDNASL.G95C
C95 is on the so-called "linker helix" of the RLC. Because of its sequence position on the helix, the spin label probably projects away from the protein. Fourteen nanoseconds of MD simulation show two modes of dynamics about different minima, each lasting ~9 and 5 ns, respectively. Given the limitations of our methods, the G95C model is most problematic because the spin label is highly exposed. In the MD simulation, the spin label appears to lie down on the surface of the protein, but this is poorly accounted for with an in vacuo model. Furthermore, such an exposed probe with few steric restrictions may have many possible conformations across a large region of conformational space with similar energies, only a few of which are revealed by 14 ns of dynamics simulation. Thus, a simulated annealing procedure was used to examine other possible minima.Simulated annealing of G95C-FDNASL
Four different conformations of the spin probe were constructed by rotating the spin probe about the cysteine sulfur and methylene carbon and locally minimizing each. These structures were then used as starting conformers for MD trajectories for 200 ps at 600 K. The three lowest-energy frames from each 600 K trajectory were chosen for subsequent energy minimization and MD simulation for 6 ns at 300 K, yielding a total of 12 trajectories. The 12 trajectories with annealed starting structures had backbone average structures unwound from their native helices. Average potential energies across the annealed trajectories are ~940 kcal/mol, whereas the original minimized structure produced a MD trajectory with an average potential energy of ~840 kcal/mol. Because the original G95C-FDNASL trajectory had the lowest energy and average conformation most congruent with the crystal structure, we used this trajectory in our analysis (Fig. 8).A105C-FDNASL
The A105C-FDNASL trajectory was computationally unstable because, apparently, the spin label and its crowded binding site were not adequately equilibrated. To remedy this, an annealing procedure was used. A 200-ps MD trajectory was simulated at 600 K, and the three lowest-energy conformations from 500 frames were chosen, called pick1, pick2, and pick3 (Fig. 9 B). All three conformations were minimized, resulting in average spin label RMS differences between the conformers of less than 1.6 Å. MD using each starting structure was simulated for 10 ns at 300 K, and each produced a convergent trajectory about its local minimum. However, upon velocity perturbation at 6 ns in the pick2 trajectory, a rotameric shift of nearby phenylalanine residue F109 into a previously unseen conformation allowed the spin label-protein system to assume a lower-energy state by nearly80
kcal/mol for the rest of the simulation. Dynamics about this new local
minimum was not only more stable, but had a calculated order parameter
of 0.78, very close to the experimentally measured order parameter of
0.8. The order parameters of the pick1 and pick3 trajectories were
~0.34 and ~0.25, respectively. This suggests that a realistic,
stable minimum was found, and we use the last 4 ns of the pick2
trajectory as the representative structure for A105C-FDNASL.
In general, we note that of all the mutant structures simulated, the
only MD trajectory that requires significant rearrangement of the side
chains for equilibration appears to be A105C, due to the "buried"
nature of the probe. In our convergence studies of C108, as in
annealing of G95C, the lowest-energy dynamics are for models that
perturb original protein structure the least. This is expected, because
we know that the attachment of a spin label preserves functionality in
the muscle fiber.
Time scales of spin label motion
Autocorrelation (Eq. 3) was performed on all trajectories over the
full 14-ns simulation, and all required multiexponential functions
(i.e., multiple correlation times) to fit each correlation function, as
expected for the persistent, quasiharmonic motions present in a MD
simulation. While correlated motions persist on many time scales, it is
useful to characterize three types of motion: a "very fast" decay
(
1 ps) reflecting the motion of the spin label atoms at 300 K; a
"fast" motion when the spin label explores the conformations about
a central dynamic mode (~50 ps), and in some cases, a "slow"
motion when the spin label samples different modes (on the order of 500 ps). To determine "very fast" and "fast" time constants, an
autocorrelation window of 50 ps was used to fit a biexponential
function. To determine "slow" time constants for bimodal
trajectories, an autocorrelation window of 1 ns was used to fit the
single-exponential function of time constant
slow. The resulting values are shown in Table
1.
The "fast" correlation times for all trajectories are near 50 ps,
and the fast-motion limit can be safely assumed if EPR spectra were to
be simulated from these MD trajectories. The "slow" rotational time
constant for D131C-FDNASL is ~210 ps. For the two modes present in
the G95C-FDNASL trajectory, if in fact they are exchanging over time,
the "slow" time constant would be even longer, on the time scale of
nanoseconds; however, as the trajectory length is limited, we cannot
say this with certainty. By comparing MD simulation with EPR studies of
free spin labels, a previous study has estimated that the viscous
effects of solvent increase c by roughly a
factor of 3 (Steinhoff and Hubbell, 1996). If this is the case, the
actual correlation time of the spin label for a bimodal trajectory may be well into the 10-ns range, a time scale in which fast motional averaging of EPR tensors cannot be assumed. This has broad implications for resolving the potentially ambiguous lineshapes attributable to time
scale and orientational disorder in an EPR spectrum. In EPR
spectroscopy, one often makes conclusions based on the assumption that
a spin label undergoes fast diffusion within a restricted cone of
motion. In some cases, bimodal exchange may be a better model, and
a detailed MD model would be an indispensable tool to help
resolve this. Also, by explicitly computing the autocorrelation directly from a trajectory, it is possible to formulate more precise simulations of EPR spectra from first principles.
Fig. 10 shows the autocorrelation of
500 ps of the D131C-FDNASL trajectory simulated at 300 and 600 K. The
autocorrelation of the first 50 ps was also done (not shown), with
t = 14.3 ps at 300 K and t = 8.0 ps at
600 K. The differences in single-exponential fits across different time
windows ( = 50 ps and = 500 ps) show the persistence
of correlated motion with increasing time scale. In both time windows,
the 600-K trajectories have rotational time constants approximately
one-half those of the 300-K trajectory, confirming the
temperature-dependent Debye relationship c ~ 
V/kT. From the asymptote of each
autocorrelation, we can discern that the 600-K trajectory, as expected,
is much more disordered. This demonstrates the difficulty in simply
increasing temperature in an attempt to allow the spin label to explore
more conformational space and thus better characterize its motion.
Whereas we can extrapolate the expected time constant of motion at a
given temperature from the Debye relation, at high temperature
the conformational space available to the probe and surrounding protein
is unrealistically large, and therefore not necessarily relevant to the
actual dynamics at ambient temperature.
|
Comparison of experimental and simulated order parameters
The order parameters calculated from both EPR spectra and MD trajectories are shown in Fig. 11. Error bars in the simulated values are the standard deviation of order parameters for 1-ns segments of the longer 14-ns trajectory. (In the case of the A105C-FDNASL, the trajectory is 10 ns long.) The order parameters of most spin-labeled mutants are excellently predicted to within 0.05, with the exception of the more exposed spin labels D131C-FDNASL and G95C-FDNASL, for which molecular dynamics trajectories predict lower order parameters than measured.
|
It is clear from the comparison that our molecular dynamics procedure provides an excellent way to perform "virtual" mutagenesis to predict spin-labeled mutants that have a high measured order parameter. Such mutants could later be expressed and used as useful experimental tools to detect the rotational motion and orientation of the proteins to which they are attached. Furthermore, by having a molecular model of the spin label dynamics at a given labeling site, it is possible to resolve EPR spectra with unprecedented detail, as the molecular dynamics trajectory contains quantitative information about orientation, amplitude, and time scale of motion, whereas most SDSL studies are necessarily qualitative.
Why do G95C-FDNASL and D131C-FDNASL predict lower order parameters? The first, most obvious, reason for this is the way order parameter is calculated from our trajectories. All order parameters are calculated with respect to the most populous mode, even for trajectories with multiple modes. Thus, the spin label's sampling of multiple minima contributes a large amount of disorder when calculated with respect to the major mode. However, because the sampling of these modes occurs on a slow time scale, the EPR spectrum is less sensitive to the motional averaging across the multiple modes, yielding a higher "apparent" order parameter. For example, consider the two dynamical modes present in the G95C-FDNASL trajectory. Each mode is stable on a time scale of several nanoseconds or perhaps longer, since the modes persist longer than can be completely characterized in the length of the simulation. The observed barrier crossing from one minimum to another is induced by a velocity perturbation. From this, we can assume that diffusion between these two modes in the simulation proceeds in the slow limit, and in this limiting case, the EPR spectrum should be a linear combination of the spectra for each mode. Because the length of the simulation does not adequately characterize the time spent around each observed minima, we are unable to determine the ratio of these components. However, an "apparent" order parameter measured from the EPR spectrum would report the highest-order mode, as reflected in the largest tensor splittings. Indeed, whereas the second mode in the trajectory has an order parameter of ~0.37, the calculated order parameter across the first 9 ns of the trajectory, which includes only the first mode, is 0.62, very close to the value of 0.60 measured from the experimental EPR spectrum.
In the case of D131C-FDNASL, the bimodal dynamics observed in the simulation explains the discrepancy between simulated and experimental order parameter. Each of the dynamical modes in the trajectory, when fit to multiple Gaussians and analyzed individually, have high order parameters (0.74, 0.87, 0.91, and 0.94), whereas the overall order parameter calculated with respect to the most populous mode has an order parameter of 0.28. The apparent order parameter of the experimental spectrum is 0.55. This, together with a molecular dynamics trajectory that displays bimodal fluctuations on a timescale of several hundred picoseconds, suggests that motional averaging on an intermediate time scale is occurring to produce the observed EPR spectrum for D131C-FDNASL.
Another reason for disagreement between the MD simulation and
experiment may be from the exclusion of solvent from the molecular dynamics model. Our initial attempts at simulating with explicit water
were unsuccessful because the computational limits of the large model
confined trajectory lengths to only several picoseconds. However,
inclusion of solvent may only have subtle effects. There is evidence to
suggest that the greatest contribution of solvent would be to simply
increase the rotational time constant of the probe due to viscous
damping, with other effects being much more subtle, and perhaps
negligible. Columbus et al. (2001) note that the energy barriers
separating solvent-exposed rotamers of the nitroxide ring of spin
labels attached to T4-lysozyme can be accounted for solely by the
viscous effects of solvent and not interactions affecting internal
energies. Steinhoff and Hubbell (1996) have speculated that with the
addition of solvent, the hydrophobic effect would increase, whereas at
the same time the van der Waals attraction would decrease, partially
canceling each others' effects. Future molecular dynamics studies that
incorporate an implicit solvent model will help further elucidate the
effect of solvent on probe dynamics. The good agreement of predicted
and experimental order parameters using an in vacuo model in this study
suggest that it is reasonable to neglect solvent except in cases where the spin label is thoroughly exposed, as in the G95C model. Exposed mutants generally have low order parameters, and because our original aim of this project was to find "virtual mutants" with high order parameter, we have seen that it is possible to neglect solvent without
catastrophic consequences.
Side chain interactions with the spin label
We know that because of surrounding residues, spin labels exhibit anisotropic dynamics. Yet SDSL is generally thought to report local structure and dynamics in an unbiased manner. Is it possible that the order parameter predicted by our MD simulation simply reflects the general topology of the site where the spin label binds, and therefore we could use local topology alone to predict good sites for binding? For each mutant trajectory, we plotted predicted order parameter versus the density of atoms within 5 Å of the spin label. There is poor correlation between the two quantities, suggesting that spin label dynamics is a more complicated story.
Do specific chemical residues near the spin label affect its dynamics? We calculated the average RMS distances and relative diffusion constants between spin label atom groups and the eight closest residues in each mutant model, but the results of these calculations showed no clear correlation with order parameter. It was observed that in the two trajectories with the highest predicted order parameter (C108 and A105C) F109 is the closest residue, but this evidence is circumstantial. Overall, we see that spin label dynamics are sufficiently complex to necessitate MD simulation.
Absolute orientation of the FDNASL in muscle fibers
Changes in the orientation of spin-labeled RLC of ~36°
with respect to the actin axis have been observed by detecting the orientational anisotropy of attached spin probes in the muscle fiber by
EPR spectroscopy (Baker et al., 1998). These data suggest that when a
muscle fiber is in rigor (myosin strongly bound to actin), the light
chain domain is predominantly at a single orientation relative to the
actin filament in oriented muscle fibers. This observation provides
another opportunity to compare experimental results with our computer
model of FDNASL attached at position 108 of the RLC. Using the
published coordinates of an actin filament decorated with chicken
skeletal myosin fragments (Mendelson and Morris, 1997), we oriented our
model of the light chain domain with respect to the actin filament.
This alignment oriented our structure in the muscle fiber reference
frame, allowing us to calculate , the angle of the nitroxide
z axis with respect to the actin filament. The aligned model
structure reported a value of = 42° whereas the
experimentally determined rigor value determined by EPR was = 38° ± 4°. This is remarkable agreement, considering several
assumptions were made in the alignment model (the decorated filament
model itself was from another species, the chicken gizzard light chain
was homology modeled onto the scallop light chain structure, and the
crystal structure probably does not reflect the inherent flexibility of
the light chain domain (Nyitrai et al., 2000; Ramachandran and Thomas,
1999; Roopnarine et al., 1998; Volkmann et al., 2000)).
| |
CONCLUSION |
|---|
|
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS AND DISCUSSION
CONCLUSION
REFERENCES
|
|---|
The order parameters calculated from MD trajectories of site-directed spin labels agree well with those determined by EPR spectroscopy, even though visual inspection of the structures surrounding these labeling sites does not provide enough information to predict probe mobilities. Thus, our procedure provides an effective way to perform "virtual" mutagenesis to predict spin-labeled residues that will display desired mobility. This can reduce the time required to synthesize, express, and assess the usefulness of constructed mutants.
MD simulation is useful as a tool to supplement EPR spectroscopy. Amplitudes of spin label motion, rotational time scales, angular distribution, and absolute orientation can be extracted from a dynamics trajectory and used to resolve ambiguous or indeterminate spectra, particularly those that may result from bimodal or multimodal dynamics. As the quality of future MD simulations improves, future studies will involve simulating accurate EPR spectra directly from molecular dynamics trajectories of spin probes, a prospect that offers a new level of detail to complement SDSL studies.
| |
ACKNOWLEDGMENTS |
|---|
We thank Adam Butensky-Bartlett, Claire Chisolm, Will Meland, and Greg Wilde, who made significant contributions to the development of this project. We also thank Octavian Cornea for help in preparation of the manuscript.
This work was supported by grants from the National Institutes of Health Grant AR32961-16, the Muscular Dystrophy Association, and the Minnesota Supercomputer Institute (to D.D.T.). L.E.W.L. was supported by a National Institutes of Health training grant and by a Doctoral Dissertation Fellowship from the University of Minnesota Graduate School.
| |
FOOTNOTES |
|---|
Address reprint requests to David D. Thomas, University of Minnesota, Biochemistry, Molecular Biology and Biophysics Department, 321 Church Street SE, 6-155 Jackson Hall, Minneapolis, Minnesota 55455. Tel.: 612-625-0957; Fax: 612-624-0632; E-mail: ddt{at}ddt.biochem.umn.edu.
Submitted January 27, 2002, and accepted for publication May 13, 2002.
Leslie E. W. LaConte and Vincent Voelz contributed equally to this work.
| |
REFERENCES |
|---|
|
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS AND DISCUSSION
CONCLUSION
REFERENCES
|
|---|
Biophys J, October 2002, p. 1854-1866, Vol. 83, No. 4
© 2002 by the Biophysical Society 0006-3495/02/10/1854/13 $2.00
This article has been cited by other articles:
 |
|
 |
|
  J. C. Klein, A. R. Burr, B. Svensson, D. J. Kennedy, J. Allingham, M. A. Titus, I. Rayment, and D. D. Thomas Actin-binding cleft closure in myosin II probed by site-directed spin labeling and pulsed EPR PNAS, September 2, 2008; 105(35): 12867 - 12872. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 K. V. Kepple, N. Patel, P. Salamon, and A. M. Segall Interactions between branched DNAs and peptide inhibitors of DNA repair Nucleic Acids Res., September 1, 2008; 36(16): 5319 - 5334. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 Y. E. Nesmelov, R. V. Agafonov, A. R. Burr, R. T. Weber, and D. D. Thomas Structure and Dynamics of the Force-Generating Domain of Myosin Probed by Multifrequency Electron Paramagnetic Resonance Biophys. J., July 1, 2008; 95(1): 247 - 256. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
S. C. DeSensi, D. P. Rangel, A. H. Beth, T. P. Lybrand, and E. J. Hustedt Simulation of Nitroxide Electron Paramagnetic Resonance Spectra from Brownian Trajectories and Molecular Dynamics Simulations Biophys. J., May 15, 2008; 94(10): 3798 - 3809. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
M. Sammalkorpi and T. Lazaridis Modeling a Spin-Labeled Fusion Peptide in a Membrane: Implications for the Interpretation of EPR Experiments Biophys. J., January 1, 2007; 92(1): 10 - 22. [Abstract] [Full Text] [PDF]
|
|
|
|
|
