Use of Molecular Dynamics Simulations in Analysis of Fluorescence-Detected Resonance Energy Transfer (FRET) Experiments

Matthew C. Zwier

Dr. Brent P. Krueger

Hope College, Holland MI 49423

March, 2003



Abstract

In the analysis of fluorescence-detected resonance energy transfer (FRET) data, it is typically assumed that both chromophores rapidly sample all orientations on a timescale that is faster than energy transfer. However, in some systems (e.g. membrane-bound regulatory proteins), one or both chromophores may be constrained to a particular orientation or range of orientations. Thus, use of the \(\ensuremath{\kappa^2}\xspace = 2/3\) approximation may not be valid. In these cases, the results of molecular dynamics (MD) simulations, taken along with fluorescence lifetime data, may yield more information than either calculation or experiment alone. We present results of MD simulations that indicate fluorophore motion is significantly restricted by conjugation with a protein (lysozyme, in this case). In addition, because our chromophores of interest contain biologically-atypical heteroatoms -- such as bromine or aromatic oxygen -- the refinement of necessary MD parameters is also presented.


Contents

Introduction

A variety of techniques capable of elucidating protein structure are in routine use, including X-ray crystallography, NMR, and fluorescence. While fluorescence techniques are more easily applied to complex systems such as transmembrane proteins, the utility of such techniques may be limited by several assumptions commonly used in order to make analysis tractable. For example, the use of fluorescence-detected resonance energy transfer (FRET) techniques as structural probes may be limited by assumptions regarding the freedom of motion of chromophores. Molecular dynamics (MD) simulations can be used to model the motion of such systems; combining experimental data with MD simulations may provide more information than is available from either technique alone.

Resonant Energy Transfer

The rate of resonant energy transfer between chromophores is a commonly-used indicator of protein structure. Because all other decay processes of a photoexcited donor chromophore should be similar in the absence and presence of an acceptor chromophore, the difference in fluorescence lifetime can be attributed to resonance energy transfer. The resonance energy transfer rate is in turn a measure of the distance between the chromophores and their relative orientations as is given in the Förster equation:

\begin{displaymath} k = \frac{3}{2} k_F \kappa^2 \left(\frac{r_F}{r}\right)^6 \end{displaymath}

In the above equation, \(k\) is the rate of energy transfer, \(r\) is the distance between the chromophores, \(r_F\) is the interchromophore distance at which energy transfer is 50% efficient, \(k_F\) is rate of transfer when \(r = r_F\), and \ensuremath{\kappa^2} is a term dependent on the relative orientation of the transition dipole moments of the chromophores. The parameters \(r\) and \ensuremath{\kappa^2} are illustrated in the figure below. Transfer slows with increasing distance and as the dipoles tend towards an orthogonal orientation.

\resizebox{\textwidth}{!}{\includegraphics{k2-blobs.eps}}

The \(r\) and \ensuremath{\kappa^2} terms are typically considered as averages over the fluorescence time frame. For systems in which both chromophores are free to assume all relative orientations within the fluorescence time frame, \(\left<\ensuremath{\kappa^2}\xspace \right> = 2/3\). For systems where motion is restricted, such as biomolecular systems, this approximation is of uncertain validity. Such restrictions are likely to be present in many biological systems of interest, e.g. transmembrane signaling proteins.

Molecular Dynamics as a Predictor of k

Molecular dynamics (MD) simulations can be used to model the motion of a system of interest and calculate \(r\), \ensuremath{\kappa^2}, and therefore \(k\). Molecular dynamics simulates the motion of large, many-atom systems using the laws of classical mechanics. This motion is governed by a force field, a combination of a potential energy function and the per-atom parameters necessary for accurate simulation of a system. The AMBER force field used in our MD simulations is shown below. It features harmonic bond stretching and angle bending potentials, a dihedral force, a 6-12 van der Waals interaction term, and a Coulombic electrostatic term.

\begin{eqnarray*} U &=& \sum_{\mathrm{bonds}} k_r (r - r_0)^2 \\ &+& \sum... ... &+& \sum_i \sum_{j \neq i} \frac{q_i q_j}{\epsilon_0 r_{i,j}} \end{eqnarray*}


The measured rate \(k\) from fluorescence lifetime experiments may be compared with \(k\) obtained from MD simulations on several systems to aid in the determination of unknown macromolecular structures, especially homooligimeric systems of known monomer structure, but unknown oligomer structure. We are developing this method using a test system of a small protein and two dyes.

Methods & Results

Protein: Lysozyme

Lysozyme is a small (<200 residues) protein readily available from chicken eggs or bacterial colonies. Additionally, it binds the dye eosin tightly and predictably. Below is a representation of hen egg white lysozyme. We produce this protein through expression by E. coli bacteria.

Space-Filling Model of Hen Egg White Lysozyme

Dye: Eosin

3D View of Eosin

Eosin, right, is an anionic dye that binds to a hydrophobic pocket of lysozyme, immobilizingthe dye. Eosin presents a challenge to MD systems because it contains two features not commonly found in biological systems. First, eosin contains covalently-bound bromine; the AMBER force field typically only uses bromine as a counterion. Second, the xanthene-like fused triple ring is completely aromatic, but contains oxygen (see resonance structure below). This aromatic oxygen does not exist in the commonly-used and well-trusted 1994 AMBER force field. Both bromine and aromatic oxygen must be parameterized for use in the AMBER system.

Resonance Forms of Eosin

Dye: DACM

Our second chromophore is DACM (below), a coumarin derivative containing an aromatic oxygen heterocycle similar to that found in eosin. This dye reacts with sulfhydryl groups to form a thioether bridge. Using this behavior, we will attach the dye to a cysteine residue in lysozyme. A mutated plasmid containing only one cysteine codon (aside from those cysteines already involved in disulfide bridges) will be used to express lysozyme for this step.

DACM Resonance Forms

Parameterization: C - Br Bonds

Aromatic carbon (CA) to bromine bonds were parameterized using a quantum mechanical (QM) method. The CA-Br bond length and CA-Br-CA bond angle were adjusted and the energy of the system calculated using a MP2/cc-pVTZ model chemistry. Parabolae were then fit to the resulting energy vs. displacement curves, as is shown below. The resulting equilibrium positions and force constants for bond and angle stretches were very similar to results published by Shearer and Cramer (J. Phys. Chem. B, 2002, 106, 5075 - 5085), and so simulation proceeded using these previously-published parameters.

QM Bond Scan Plot

Parameterization: Aromatic Oxygen

Structure of Fluorescein

The aromatic oxygen heterocycle present in eosin and DACM posed a greater challenge. There was no straightforward way to separate bond stretching and angle bending contributions from aromatic oxygen (OA) and aromatic carbon (CA) in a quantum-mechanical energy surface simulation. Thus, an alternative method was employed. Fluorescein, right, was used as the model dye; it contains the features of eosin without the electron-rich bromine substituents, which slow QM calculations considerably.

OA out of plane bend

A quantum mechanical calculation of the normal vibrational modes of the system was performed, normal modes strongly dependant on OA movement were identified, and OA parameters adjusted until the normal mode frequencies produced by AMBER were sufficiently close to those produced by the quantum mechanical model. Van der Waal’s parameters for OA were assumed to be similar to ether-type oxygen (OS). Bond, angle, and dihedral parameters were optimized with the above method. Special emphasis was placed on the dihedral parameters under the assumption that a disproportionate out-of-plane bend of the aromatic oxygen in eosin or DACM would produce exaggerated fluctuations in the dipole moments of the chromophores, which would in turn lead to inaccurate results for MD-calculated \ensuremath{\kappa^2} values. One such dihedral parameter, the out-of-plane OA bend, is illustrated at the left.

Optimized OA Parameters

The following table demonstrates the optimized OA parameters:

Original AMBER Optimized AMBER

Description

QM Freq. Frequency Ratio Frequency Ratio
Out-of-plane ring twist 115.5 139.9 1.21 140.2 1.21

Out-of-plane OA bend

277.6 221.4 0.80 280.4 1.01

In-plane ring bend

697.3 696.8 1.00 704.2 1.01

Aromatic ring C stretch

1681.3 1790.9 1.07 1790.9 1.07
All frequencies are given as wavenumbers ( \ensuremath{\mathrm{cm^{-1}}}). QM frequency refers to the normal mode frequency as calculated by Gaussian 98 with the B3LYP/6-31G* model chemistry. Original AMBER frequency refers to the normal mode frequency as calculated using ether-type oxygen parameters, and optimized AMBER frequency refers to the normal mode frequency calculated with adjusted OA parameters. Ratios are (AMBER freq)/(QM freq).

Molecular Dynamics

A molecular dynamics simulation was then performed on a system containing lysozyme, eosin, and DACM (below) solvated in water in a periodic box with sides of length 12 Å. Electrostatic interactions were treated with particle-mesh Ewald (PME). The production run was preceded by 500 steps of conjugate-gradient minimization, 2500 steps of constant-volume, constant-temperature equilibration, and 5000 steps of constant-pressure, constant-temperature equilibration. A trajectory of 925 ps was analyzed. Click on the image to see a ~ 10 MB movie of the trajectory.

Image chx+dacm+eosin-lores

Results: \ensuremath{\kappa^2}

Simple trigonometric treatment of the dipole moment vectors of the chromophores yields a value for \ensuremath{\kappa^2} at each time step, plotted below. The average value of \ensuremath{\kappa^2} is 0.56, or 16% lower than the \(\left<\ensuremath{\kappa^2}\xspace \right> = 2/3\) approximation.

Kappa Squared vs. Time

Results: k

From the values of \ensuremath{\kappa^2} calculated above, a relative predicted energy transfer rate was calculated as \(\ensuremath{\kappa^2}\xspace /r^6\), which is shown in blue in the plot below. The corresponding curve with \ensuremath{\kappa^2} fixed at 2/3 is shown in red. The relative energy transfer rate using the calculated \ensuremath{\kappa^2} is predicted to be 29% lower than that expected by the \(\left<\ensuremath{\kappa^2}\xspace \right> = 2/3\) approximation.

Energy Transfer Rate vs. Time

Conclusions & Discussion

We have demonstrated a feasible molecular dynamics-based method for calculating resonant energy transfer rates in systems of known structure on a system consisting of lysozyme, eosin, and DACM. AMBER force field parameters for the biologically-atypical atoms (bromine and aromatic oxygen) in our chromophores were verified and developed, respectively.

In the lysozyme-eosin-DACM system, eosin is immobile, and one would expect deviation from the \(\left<\ensuremath{\kappa^2}\xspace \right> = 2/3\) approximation. In fact, in a short molecular dynamics run, \(\left<\ensuremath{\kappa^2}\xspace \right>\) was calculated as 0.56, and the corresponding energy transfer rate is ~ 30% lower than that predicted by the \(\left<\ensuremath{\kappa^2}\xspace \right> = 2/3\) approximation.

In the future, we will measure the resonant energy transfer rate of the lysozyme-eosin-DACM system using fluorescence lifetime spectroscopy. This will allow comparison to and refinement of the computational model. Eventually, the model will be applied to transmembrane signaling proteins, in which one can expect significant deviation from the \(\left<\ensuremath{\kappa^2}\xspace \right> = 2/3\) approximation.

Acknowledgements

Quantum mechanical calculations were performed using Gaussian 98, and molecular dynamics simulations were carried out using AMBER 7. Protein system illustrations were generated with VMD (Visual Molecular Dymanics) from the Schulten group at the University of Illinois at Urbana-Champaign. Small three-dimensional molecular illustrations were generated using POVChem. Both protein and small molecular system images were rendered with POV-Ray.

Many calculations were carried out at the National Center for Supercomputing Applications (NCSA), which we thank for providing supercomputer time. Along with the Hope College Department of Chemistry, we also thank the National Science Foundation through the Research Experiences for Undergraduates program and Research Corporation for funding.

We extend our gratitude to Christopher Cramer, Junmei Wang, and Xanthipe Jordanides for helpful correspondence.