Article Information

• PDF (713 kb)
• Full Text with Large Figures
• Cited by in Scopus (14)

PubMed

• Articles by Gunnar F. Schröder
• Articles by Helmut Grubmüller

Related Articles

• Correction

  Correction Biophysical Journal, Volume 87, Issue 2, 1 August 2004, Pages 1398 Full Text | PDF (43 kb)

• Total Internal Reflection Fluorescence Correlation Spectrosc...

  Total Internal Reflection Fluorescence Correlation Spectroscopy: Effects of Lateral Diffusion and Surface-Generated Fluorescence Biophysical Journal, Volume 95, Issue 1, 1 July 2008, Pages 390-399 Jonas Ries, Eugene P. Petrov and Petra Schwille Abstract Fluorescence correlation spectroscopy with total internal reflection excitation (TIR-FCS) is a promising method with emerging biological applications for measuring binding dynamics of fluorescent molecules to a planar substrate as well as diffusion coefficients and concentrations at the interface. Models for correlation functions proposed so far are rather approximate for most conditions, since they neglect lateral diffusion of fluorophores. Here we propose accurate extensions of previously published models for axial correlation functions taking into account lateral diffusion through detection profiles realized in typical experiments. In addition, we consider the effects of surface-generated emission in objective-based TIR-FCS. The expressions for correlation functions presented here will facilitate quantitative and accurate measurements with TIR-FCS. Abstract | Full Text | PDF (517 kb)

• Pattern Formation by Bacteria-Driven Flow

  Pattern Formation by Bacteria-Driven Flow Biophysical Journal, Volume 88, Issue 4, 1 April 2005, Pages 2525-2529 N.G. Cogan and Charles W. Wolgemuth Abstract Some marine bacterial species form mucosal layers, called veils, on sulfidic marine sediment. The bacteria attached to this veil actively swim and exert force on the surrounding fluid. The bacteria can break free of the veil and swim, chemotacting back to the veil. Over time the veil forms holes arranged in a hexagonal pattern. Motivated by this system, we present a simplified model to describe pattern formation induced by force-generating bodies embedded in a layer surrounded by fluid. When the bacteria break free of the layer, they are advected by the flow and diffuse. Competition between the fluid flow generated by the embedded bacteria and diffusion of the swimmers leads to a novel instability that drives bacterial aggregation. Analytic and numeric analysis of this system correctly defines the length scale and developmental timescale for the biological system. Similar flow dynamics may also play a role in other biological systems such as encrusting bryozoan colonies. Abstract | Full Text | PDF (141 kb)

• …more

Copyright © 2005 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 89, Issue 6, 3757-3770, 1 December 2005

doi:10.1529/biophysj.105.069500

Biophysical Theory and Modeling

Simulation of Fluorescence Anisotropy Experiments: Probing Protein Dynamics

Gunnar F. Schröder^*, Ulrike Alexiev^† and Helmut Grubmüller^*, ,  

^* Theoretical and Computational Biophysics Department, Max-Planck-Institute for Biophysical Chemistry, Göttingen, Germany
^† Department of Physics, Freie Universität Berlin, Berlin, Germany

Address reprint requests to Helmut Grubmüller.

Protein motions, particularly their conformational dynamics, regulate and often constitute protein function. Therefore, a large variety of experimental and theoretical techniques aim at probing internal dynamics of proteins, with a particular focus at the picosecond to microsecond timescale. Nuclear magnetic resonance, electron paramagnetic resonance, neutron scattering, as well as fluorescence depolarization experiments have indeed provided much insights in this respect.

Fluorescence spectroscopy, in particular, in combination with site-directed fluorescent labeling has become an established tool to investigate the dynamics and interactions of biomolecules. Fluorescence resonance energy transfer experiments allow us to determine intramolecular distances in biomolecules 1,2, recently even with millisecond time resolution 3,4,5. Time-resolved fluorescence anisotropy experiments provide information on both, mobility and dynamics of a fluorophore 6,7. Because this probed dynamics of the dye is affected by the motion of the protein fragment to which it is attached, the anisotropy decay yields information on the protein structure and conformational changes as well as on the protein flexibility. Recently, this approach has been successfully applied to surface loop dynamics of rhodopsin and bacteriorhodopsin (bR), where an anisotropy decay component could be assigned to the loop motion, which allowed for the identification of protein conformational changes 8.

To extract structural information from the anisotropy decay curves is not straightforward, however. Several empirical models have therefore been proposed to facilitate the interpretation 9,10,11,12. Accordingly, many results obtained by these experiments depend on the particular choice of model used. Furthermore, each model rests on certain assumptions that are not always easy to verify. A first-principles simulation approach, in contrast, would enable one to drop these assumptions and thus to provide more accurate interpretations of the experiments at the molecular level. To this aim, we have carried out molecular dynamics (MD) simulations of the complete experimental system including the protein, the protein-attached dye and an explicit solvent environment. These simulations should allow us to extract the individual contributions to the depolarization and to analyze the dye-protein interactions at atomic detail.

In a similar spirit, but at faster timescales, MD simulations of tyrosine- and tryptophan-containing proteins have been used to predict the fluorescence anisotropy decay function 13,14,15,16,17. Furthermore, anisotropy decay measurements of a dye attached to flexible protein regions have been correlated with side-chain dynamics from MD simulations 18. Also, simulations of free fluorophores in a solvent reproduced temperature and solvent dependence of the experimental fluorescency anisotropy 19,20. In particular, the rotational diffusion of tryptophan in water has been simulated and its dependency on different water models discussed 21,22.

Here, and in contrast to the previous studies, the focus is on the interaction and dynamic coupling between a particular protein fragment and the attached dye. In addition, determination of dye conformation often is a key to the interpretation of fluorescence spectroscopy experiments. We will show that dye conformations can also be obtained from MD simulations, which therefore provide valuable information complementary to the experiment.

The system studied here is a fragment of the heptahelical transmembrane protein bR, which shares the heptahelical bundle motif with the large class of G-protein coupled receptors (GPCRs) and is considered a prototype for GPCRs. In general, GPCRs participate within many cell signaling pathways by conversion of external stimuli into intracellular responses. The dynamics and flexibilities of surface-exposed protein loops of these GPCRs have been shown to play a crucial role in activation and molecular recognition of the receptor, as G-proteins transduce the conformational change of the activated receptor into a regulated effector function (for recent reviews, see Sakmar et al. 23 and Marshall 24). Investigation of the local protein surface dynamics and flexibility can thus provide much insight into the mechanism of receptor activation and signal protein binding to the protein surface. Recent picosecond time-resolved fluorescence anisotropy decay measurements 8 have addressed the dynamics of these loops.

Here we report on a combined MD simulation and experimental approach, where both, measurements and simulations have been devised to match each other as closely as possible. In the experiments, and in the simulations we used the C₅ maleimide of the fluorescent dye Alexa488 covalently bound to a single cysteine residue introduced at position 35 of bR to probe the flexibility of the loop connecting the first two transmembrane helices A and B. In close analogy to the experimental design used for the assignment of the loop component 8, the AB fragment of bR has been studied at first to characterize the intrinsic dynamics of the AB loop. In a next step, comparison with loop dynamics of the intact protein is feasible and is expected to provide information on the interactions between the loops.

Our aim in this study is to gain insight into fluorescence anisotropy experiments by providing an interpretation of the experiment in atomic detail. We particularly address the following questions: Which processes influence the reorientational dynamics of the dye? What are the individual contributions to the observed depolarization? How does one extract information on the protein conformational dynamics from the anisotropy decay curve? Vice versa our simulations should enable us to study to what extent does the attached dye affect the unperturbed loop dynamics. This effect is commonly—and necessarily—assumed to be small. This study offers the chance to test this assumption.

Here we summarize the basic equations that will be used to calculate anisotropy decay curves from molecular dynamics simulations as well as from experiments. In the experiments, dye molecules are excited by a short polarized laser pulse. After a period of time t, during which the dye has undergone rotational diffusion, the polarization of the reemitted light is measured. The fluorescence anisotropy r(t) at time t after excitation of the fluorophore is defined as

(1)

(2)

For the case at hand, where the dye has been attached to a protein, the motional freedom of the dye is restricted. A common model to describe such restricted rotational diffusion, which we will also use as a reference, is the “wobbling-in-a-cone” model 9. In this model, the transition dipole is assumed to diffuse freely inside a cone (Fig. 1, top) fixed within the molecular frame, and the anisotropy r(t) can be approximated by

(3)

Display large version of this figure

Figure 1

(A) Wobbling-in-a-cone model. The transition dipole moment μ(t) is assumed to freely diffuse inside a cone with the semicone angle θ_max. (B) Cone-in-a-cone model. The two top figures indicate a motion, which is due to the protein flexibility, superimposed to the wobbling-in-a-cone model. The bottom figure visualizes how the protein flexibility is described in the cone-in-a-cone model by the protein cone, within which the dye cone freely diffuses.

Here, r₀=0.4 P₂(cos λ), and λ is the angle between the absorption and emission dipole moment fixed in the coordinate frame of the dye. A_∞ is a parameter describing the degree of motional restriction and is related to the semicone angle θ_max (Figure 1A) by

(4)

Note that a small cone angle gives rise to a large value of A_∞; vice versa, A_∞ vanishes for isotropic rotational diffusion. Additionally considering an isotropic overall tumbling motion of the dye-protein complex as a whole, described by a rotational correlation time ϕ_G, the anisotropy of the protein-attached dye is given by

(5)

The orientational dynamics of the dye is also affected by the local flexibility of the protein, which is of particular interest here. The resulting surface changes affect cone orientation and possibly, also cone angle (Figure 1B, top). In a simple model 14, also these protein-induced changes are described by a restricted motion of the dye cone within a (second) cone (Figure 1B, bottom). Assuming that the local dye motion and the protein dynamics are uncoupled, this effect gives rise to a second decay factor,

(6)

It is therefore referred to here as the “cone-in-a-cone” model.

The preparation and expression of the bR mutant S35C in Halobacterium salinarium, in which cysteine replaces Ser-35, has been reported 27. Solubilization of bR membrane fragments and regeneration of bR in DMPC/CHAPS-micelles were performed as reported 27. The regeneration procedure was modified according to protocols described 28. Briefly, bR membrane fragments were solubilized in 1.2% (w/v) SDS and 1mM DTT. After the retinal band is shifted completely to ∼380nm the excess SDS was removed by chromatography on Sephadex-G25, preequilibrated, and eluted with 0.1% SDS in 5mM sodium phosphate buffer at pH 6. The regeneration was started by adding 2% CHAPS/2% DMPC, 20mM sodium phosphate buffer, pH 6.2, and 400mM KCl. After regeneration excess lipids, detergents and buffer were removed by chromatography on Sephadex-G25, preequilibrated, and eluted with 0.1% CHAPS and 150mM KCl. The retinal band is shifted back completely to 550nm, resulting in 96–100% regeneration, using an estimated extinction coefficient of ϵ₅₅₀≈52,000M⁻¹cm⁻¹. Labeling of bR with Alexa488 and the determination of the labeling stoichiometry were performed as described 29 for 5- (iodoacetamido)-fluorescein and adapted for Alexa488C₅ maleimide (Molecular Probes, Eugene, OR). The reaction of 50μM bR with 2mM fluorescent label was carried out at room temperature in 150mM KCl, 100mM Tris buffer, pH 8.0, 30mM EDTA, 200μM DTT. Excess reagents were removed by chromatography on Sephadex G-25 in 0.1 A 10-mM stock solution of the respective label was used. The labeling stoichiometry was calculated using

(7)

Display large version of this figure

Figure 2

(A) SDS-PAGE of fluorescently labeled bacteriorhodopsin S35C-Alexa488. (Lanes 1 and 4) Molecular weight marker, (lane 2) undigested bacteriorhodopsin S35C-Alexa488, (lane 3) partially digested bacteriorhodopsin S35C-Alexa488. The AB fragment in lane 3 is marked by a white arrow. (B) Secondary structural model of the AB-helix fragment of bacteriorhodopsin. The amino acids are shown in single letter code.

Display large version of this figure

Figure 3

Absorption (black) and emission (green) spectra of the fluorescent dye Alexa488 covalently bound to a single cysteine on the surface of bacteriorhodopsin in position 35. The absorption spectrum unlabeled bacteriorhodopsin (pink) is shown for comparison. The spectra are scaled. The absorbance and emission maxima of Alexa488 are marked. The labeling stoichiometry Alexa/bR was 0.63. Conditions are: 0.1% CHAPS, 150mM KCl for the absorbance spectra; methanol for the fluorescence emission spectrum of the fragment (1–72) S35C-Alexa488. The inset shows the chemical structure of Alexa488; the arrow indicates the transition dipole moment of the chromophore.

The proteolytically cleaved fragments of the fluorescently labeled bR mutant S35C containing residues 1–72 were used as a model 2-helix system (Figure 2B) to analyze the anisotropic behavior of the fluorescent dye in its protein bound state 8. Before time-resolved fluorescence measurements, the AB-helix fragment was concentrated and dissolved in methanol.

The fluorescence anisotropy decays were measured employing a tunable Ti:sapphire laser/microchannel plate-based single-photon counting apparatus with picosecond time resolution. A detailed account of the experimental setup has already been provided 8. This method allows to measure the diffusional dynamics of the loops directly on the pico- to nanosecond timescale. The fluorescence decay profiles, I_‖(t) and I_⊥(t) or I_‖(t)+2I_⊥(t), and the time-resolved anisotropy as given by Eq. (1) were analyzed using the software package Global Unlimited V2.2 (Laboratory for Fluorescence Dynamics, University of Illinois, Urbana, IL). The time course of the fluorescence was deconvoluted and fitted with a sum of exponentials

(8)

The anisotropy decay was fitted with a model function comprising a sum of three exponentials as well as with Eq. (6).

Three simulation systems were constructed. The first two systems contained one Alexa488 dye molecule (see inset of Fig. 3) solvated in methanol and water, respectively. The third system comprised the AB fragment of bR with the Alexa488 dye attached to position S35C, solvated in methanol. All MD simulations were performed using the GROMACS simulation software 31 with the united-atom GROMACS force field, which describes nonpolar hydrogens implicitly via compound atoms. The compound atom approximation could potentially underestimate C-H/aromatic interactions between the protein and the dye. However, we do not expect this effect to significantly influence our results. The SPC 32 water model was used. The methanol parameters were taken from the GROMACS force field. The force-field parameters for the dye have been determined as described below. All systems were energy minimized to obtain the starting configuration for the simulations. Periodic boundary conditions were applied. In the free-dye simulations the solvent and the solute were jointly coupled to an external temperature bath of 300K with a relaxation time of 0.1ps 33. For the protein-dye system the solvent was separately coupled to a heat bath with the same parameters. In all simulations the system was weakly coupled to a pressure bath of 1atm with isotropic scaling and a relaxation time constant τ_p=1ps. Bond lengths were constrained to their equilibrium lengths using the LINCS algorithm 34. This allowed for a 2-fs time step using the leapfrog integration scheme. For the Lennard-Jones interactions, a cutoff distance of 1.0nm was applied. Electrostatic interactions between charge groups at a distance <1nm were calculated explicitly, and the long-range electrostatic interactions were calculated using the particle-mesh Ewald method 35 with a grid spacing of 0.12nm and a fourth-order spline interpolation. Coordinates of all atoms were saved every 1ps for further analysis. The initial structure of the AB fragment (residues 8–71) of bacteriorhodopsin was taken from the crystal structure, Protein Data Bank entry 1AP9 36. It has been shown that the AB fragment solvated in an organic solvent (methanol/chloroform (1:1)) adopts a conformation similar to its structure in bR 37.

Because the motional restriction of the dye due to the protein surface is mainly determined by steric hindrances and electrostatic interactions, we paid particular attention to those force-field parameters that sensitively affect these quantities, i.e., van der Waals parameters and partial charges. The van der Waals parameters are relatively insensitive to the chemical environment and were thus taken from corresponding atom types from the GROMACS force field. For the aliphatic linker region, parameterization was straightforward because this region is chemically sufficiently similar to aliphatic groups in the GROMACS force field. Similarly, the force constants describing the chemical bonds within the headgroup of the dye have been adapted from the GROMACS force field. Due to the extended π-electron system of the dye, this simple estimate is not expected to yield accurate force constants; however, as the internal vibrations of the dye are largely uncoupled from the overall diffusive dye motion of interest, this inaccuracy is not expected to affect our results. For the free dye, this assumption has been verified by test simulations with varying force parameters (data not shown).

The partial charges for the dye headgroup are critical and therefore have been determined from QM calculations. In all QM calculations described below, only the headgroup of the dye is considered, with the linker replaced by a methyl group. Here, the calculation is complicated by the fact that, to simulate the experiment properly, the molecular dynamics between absorption and emission have to be studied. Therefore, the dye parameters have to be determined for the first excited state.

To that aim, we followed a hybrid strategy. In a first step, the ground state of this molecule was calculated using density functional theory (DFT) as implemented in the DMol program 38 with the DNP basis set and the BLYP functional. Then both the ground and first excited states were calculated using CIS/STO-3G with the GAUSSIAN program (ab initio) 39 and also, as a check, at the semiempirical level PM3 with the program MOPAC 40. Atomic charges from all calculations were obtained by fitting to the electrostatic potential (“ESP charges”) 41. The differences of the charges between the ground and first excited state obtained from the ab initio and semiempirical calculations were small; we therefore took their mean values to minimize the complementary systematic errors inherent to both methods. These differences were then added to the ground-state DFT charges. To assure compatibility to the GROMACS protein force field, partial charges for all 20 amino acids were calculated using the same DFT level. From comparison to the respective GROMACS partial charges, a mean scaling factor of 0.7 was obtained and applied to the dye charges, including a shift to retain the correct total charge of −2e. Finally, the slightly unsymmetric charges were symmetrized, according to the symmetry of the dye, to resemble the ensemble average in a mean field approach. Force-field parameters of the bond lengths and angles were obtained from the ab initio geometry calculation.

The Alexa488 (Molecular Probes) dye is only available as a mixture of 5- and 6-isomers (cf. inset of Fig. 3); we restricted our simulations to the 5-isomer, assuming only minor dependence of the dye dynamics on the choice of the isomer. The parameter files for this dye can be downloaded from our website (http://www.mpibpc.gwdg.de/groups/grubmueller/probflex/).

For efficient sampling of the conformational space of the dye, simulations at 1000K with implicit solvent were carried out. The central oxygen atom in the headgroup of the dye (cf. inset of Fig. 3) was chosen to represent the conformation of the dye, and the probability distribution p_high of this atom at high temperature was obtained from the simulation. Because the equilibrium distribution of the dye is predominantly governed by the electrostatic interactions between the dye and the protein, this interaction was calculated for each conformation visited in the simulation. To account for the dielectric properties of the thin methanol layer between the dye and the protein, we derived a corrected dielectric coefficient by scaling a distance-dependent dielectric coefficient obtained for water 42. Following this approach, and using the distance between the dye and the protein of ∼0.5–0.8nm as observed in the simulation, a dielectric coefficient of ϵ=10 was obtained and used for the calculation of the electrostatic energy. This energy E_Coulomb served to estimate the probability distribution p at room temperature from the high temperature distribution p_high via the Boltzmann factor,

(9)

with T_high=1200K and T=300K. We note that the obtained probability distribution p (Figure 7C), is expected to provide only a rough estimate. Certainly, solvent effects will not fully be described by the implicit solvent treatment, which was necessary to obtain sufficient sampling.

The correlation between the dye and the protein was calculated using a modified version of the LMLA (localized mean by linear assignment) algorithm 43. This algorithm serves to obtain curvilinear principal coordinates from molecular ensembles of conformations, represented as sets of points x⁽ⁱ⁾ in the 3N-dimensional configurational space of an N-atomic molecule. The main idea of this approach is that a large number of randomly chosen k-tuples of structures probe the shape the molecular ensemble and, therefore, by properly averaging over all these k-tuples, information on the “average shape” of the whole ensemble can be obtained. This information is represented in terms of k “prototypic structures” (PS) in the conformational space, which are positioned along the largest extension of the molecular ensemble, thus capturing the main conformational changes of the system.

The PS a_j are calculated as localized ensemble averages

(10)

(11)

Here we are not interested in a description of the conformational changes of the whole system, i.e., dye and protein. Rather, we wish to analyze which motions of the protein are probed by the motion of the dye attached to it, i.e., which collective mode of motion is best correlated with the configurational subspace defined by the dye. To this aim, the LMLA algorithm is modified such that the linear assignment of Eq. (11) is done only within the subspace of the dye coordinates,

(12)

(13)

Here, k=3 was found to describe the nonlinear motions of interest sufficiently accurately; therefore, three subspace-determined PS (conformations of the dye-protein system) were calculated from the ensemble generated by the explicit solvent MD simulation at 300K.

To quantify the contribution of each single protein atom to the obtained correlation with the dye, we calculated the root mean square fluctuation (rmsf) σ_i of the obtained for each atom i,

To characterize the fluorescence anisotropy decay on different timescales, we calculated a position-dependent contribution to the depolarization of the dye. To this aim, only five degrees of freedom of the dye trajectory were considered: The position of the dye, represented by the center of mass x_m of the headgroup of the dye, and the normalized vector of the transition dipole moment μ(t). To gain information on the timescale of the respective correlations, we compared the original MD trajectory to a smoothed trajectory, where the fast fluctuations were filtered out and only the slow components of the dye dynamics remained. Only the transition dipole vector was low-pass filtered (yielding μ_s(t)) using a Gaussian kernel with a mean±SD σ=40ps, whereas the x_m trajectory remained unchanged.

To obtain a spatially resolved picture, the contributions ξ(t′)=P₂[μ(t′)·μ(t′+t)] to the anisotropy r(t) at t=50ps (cf. Eq. (2)) were calculated. The position-dependent contribution to the depolarization , which is a measure for the orientational flexibility at a given position x of the dye, is then obtained by

(15)

The orientation distribution (cf. Fig. 14), represented as a histogram on the surface of a sphere, has been built up of 500 cones pointing away from the average center-of-mass position of the headgroup of the dye. Their orientation is defined by 500 randomly chosen unit vectors. Their lengths n_j represent the number of transition dipole orientations that fall into the respective directions,

(16)

To assess the orientational distribution of the dye in conformation A within the loop frame, all structures from the trajectory were aligned to minimize the root mean square deviation rmsd of the loop residues 30–42. For the distribution in conformation B, all structures were aligned to minimize the rmsd of the helical residues (10– 29, 43–61).

The anisotropy decay r(t) due to the change in the orientation of the transition dipole moment of the dye between absorption and emission has been calculated according to Eq. (2). To this aim, the time-dependent absorption and emission dipole moments, μ_a(t)=M(t)·μ_a and μ_e(t)=M(t)·μ_e, had to be derived from the MD simulations. Here, μ_a and μ_e are the dipole moment vectors in the coordinate frame of the moving dye and M(t) is a rotation matrix that transform the dye coordinate frame into the lab frame. The absorption dipole moment μ_a in the dye frame has been obtained from the CIS calculation of the dye described above and is oriented along the three-ring system of the chromophore as shown in the inset of Fig. 3 (arrow). The rotation matrix M(t) is calculated for each snapshot of the MD trajectory (in 1-ps steps) from the instantaneous structure of the dye, thus yielding a trajectory of the absorption dipole vector μ_a(t). The emission dipole moment in the dye frame, μ_e, was not calculated explicitly. Rather, the angles λ=10° and λ=13° between μ_a and μ_e were determined from the measured initial anisotropies r₀=0.37=0.4P₂(cos λ) and r₀=0.34 for the label free in solution and for the label-loop construct, respectively, and found to be small. In this case, λ enters into r(t) as an overall scaling factor, and thus the normalized decay curve studied here is not affected. Therefore, μ_e(t) was replaced by μ_a(t).

The time average used to calculate the anisotropy decay r(t) (cf. Eq. (2)), obtained from one trajectory, implies a statistical error. To estimate this error, Brownian dynamics simulations have been carried out for the cone-in-a-cone model. In these simulations, 230 rotational diffusion trajectories of a single normalized vector, representing the transition dipole moment of the dye within the dye cone, as well as the diffusion of the dye-cone axis within the protein cone have been calculated. The lengths of these trajectories (16ns) have been chosen identical to the length of the MD trajectory of the full dye-protein system. The rotational diffusion coefficients and semicone angles for the dye and the protein cone (0.0018ps⁻¹, 0.0003ps⁻¹ and 45°, 50°, respectively) were chosen such as to match those extracted from the MD simulation (Table 1). From each of these 230 trajectories r(t) has been calculated as described above, and for each r(t), cone-in-a-cone parameters were obtained. From their variances, the statistical error of the cone-in-a-cone parameters derived from the MD simulation was estimated (Table 1).

In a first step, we studied if and to what extent MD simulations allow us to predict the fluorescence anisotropy of a dye and if the used dye and solvent force-field parameters are accurate enough for this purpose. To this aim, simulations of the free dye in methanol and water were carried out and the fluorescence anisotropy calculated from the simulations was compared to experiment. The Alexa488 dye was used both in the experiment described further below, as well as in the simulations.

One Alexa488 molecule was simulated in a box of 4673 water molecules for 8ns and in a box of 1317 methanol molecules for 13ns. The dye was free to undergo translational and rotational diffusion within the (periodic) simulation volume. Figure 4A shows the calculated anisotropy as solid lines (green, methanol; blue, water). Fits to a single exponential function yield a rotational correlation time of ϕ =51ps for the dye in water and ϕ=86ps in methanol. Figure 4BC, show the measured time-resolved fluorescence intensities of Alexa488 in aqueous solution and in methanol at 300K. The corresponding anisotropy decay curves, calculated from the two components of fluorescence decay I_‖(t) and I_⊥(t) according to Eq. (1), are shown as dashed lines in Figure 4A. The fluorescence lifetime in water was fitted with two exponentials, τ₁=1.97±0.03ns and τ₂=8.68±0.03ns and the anisotropy decay was fitted with one exponential, yielding a rotational correlation time ϕ=170ps. The fluorescence lifetime and the anisotropy decay of Alexa488 in methanol were fitted with single exponentials, τ=3.82±0.02ns and ϕ=210ps, respectively. The residuals of the fit for I_‖(t) and I_⊥(t) are shown as insets in Figure 4AB, respectively. The final χ² of the fit was 1.02 for Alexa488 in water and 1.007 in methanol. IRF is the instrument response function with a full width at half maximum of 48ps.

Display large version of this figure

Figure 4

(A) Calculated fluorescence anisotropy decays (solid lines) for Alexa488 in water (blue) and methanol (green) are compared to the measured anisotropy decays (dashed lines). (B and C) Time courses of the fluorescence intensities I_‖(t) (red) and I_⊥(t) (blue) in water (B) and methanol (C). IRF is the instrument response function. The residuals of the fit for I_‖(t) and I_⊥(t) are shown as insets. The excitation was at 470nm. The fluorescence emission was detected after passing through a cutoff color glass filter OG 515.

As can be seen, the correlation times in the simulation are smaller than in the experiment by about a factor of 3. Further studies 43 suggested this effect to be due to the solvent force fields. In particular, self-diffusion coefficients and dielectric relaxation times of solvent molecules obtained from MD simulations were compared to experiment and found to deviate systematically from each other. This deviation could explain the too-fast rotational diffusion of the dye molecule in the simulation. Nevertheless, the details of this effect are not fully understood yet and need to be further investigated.

Both, in the simulation as well as in the experiment the dye shows faster rotational diffusion in water than in methanol. This behavior was unexpected because the viscosity of water (1.002×10⁻³Pa s (at 293K)) is larger than that of methanol (0.587×10⁻³Pa s (at 293K)). Furthermore, a similar experiment with fluorescein dyes shows the normal behavior; with measured rotational correlation times of 140ps in methanol and 170ps in water (data not shown).

To explain this inverse solvent effect for Alexa488, the structure of the dye in the simulation was analyzed in more detail. Visual inspection suggested that the extension of the dye strongly depends on the solvent. Figure 5C quantifies this extension by the distance d defined in Figure 5A, for both simulations in methanol (green curve) and water (blue curve). As can be seen, d fluctuates between 0.7 and 1.6nm. The representative dye conformations are shown in Figure 5AB, for d=1.6nm and d=0.7nm, respectively. As the structures show, the headgroup of the dye is rather stiff, thus the change of the length is only due to the flexible linker, a hydrophobic chain. The average length of the dye in methanol of ∼1.5nm is clearly larger than that in water of ∼1.0nm (thick horizontal lines).

Display large version of this figure

Figure 5

Extension of the Alexa488 dye molecule solvated in water and in methanol, as quantified by the distance d defined in panel A (arrow). This distance is plotted versus time for the water (blue) and methanol (green) simulations. (B) Dominant conformation of the dye solvated in water, with a distance d=0.8nm.

This observation can be explained in terms of the water-exposed hydrophobic surface. In water, the hydrophobic chain minimizes the water-accessible surface by coiling up and thereby reducing the effective size of the molecule. In the more hydrophobic solvent, the hydrophobic forces are absent, and the chain remains in the entropically more favorable extended conformation. We suggest this as an explanation of the observed faster rotational diffusion of Alexa488 in water than in methanol. Furthermore, the lack of this flexible hydrophobic chain in fluorescein explains, why in this case this inverse solvent effect was not observed. These results show that the simulations of the dye are capable of describing the subtle difference in the rotational diffusion of the dye solved in water and methanol.

Having passed these initial tests, we simulated the Alexa488 dye attached to the S35C position of the AB-helix fragment (residues 8–71) of bR solvated in 18,752 methanol molecules and including two sodium ions for 26ns. Fig. 6 shows the rmsd of the backbone atoms of the helical part (residues 10–29 and 43–61). The rmsd reaches a relatively low mean rmsd value of 0.12nm after only ∼20ps, which indicates that the α-helical structure remains very stable during the simulation. Furthermore, the root mean square fluctuation (rmsf) around the mean structure of ∼0.02nm indicates a low flexibility of the protein. In contrast, the loop connecting the two helices shows a much higher flexibility than the helical regions, which will be discussed further below.

Display large version of this figure

Figure 6

Root mean square deviation (rmsd) of the backbone atoms of the helical residues (residues 10–29 and 43–61) from the initial x-ray structure during the first 10ns of the simulation.

The initial dye geometry was chosen in extended conformation, detached and pointing away from the surface. After ∼500ps, the dye reached conformation A, shown in Figure 7A, and remained loosely bound via noncovalent interactions to the protein surface. It then detached from the surface and flipped back toward the other side into conformation B, shown in Figure 7B, where it remained for the rest of the simulation. In conformation B, the dye is much less flexible than in A, as quantified by the rmsf of 0.15 and 0.23nm, respectively. In conformation A, the dye adopts two conformational substates, shown in Fig. 8 (“up” and “down”). During the residence time of 16ns in conformation A, the dye flipped back and forth three times between these two substates, which significantly contributes to the mobility of the dye. This complex and hierarchical motion is apparently insufficiently described by a simple wobbling-in-a-cone model.

Display large version of this figure

Figure 7

Snapshots of the bR protein fragment with covalently bound dye in different conformations (side view), as obtained from the MD simulation. During the first 16ns, the dye adopts conformation A (A) and is highly mobile. For the rest of the simulation, the dye remains in conformation B (B) and is much less flexible, as it is more tightly attached to the two helices. (C) The probability distribution of the dye calculated from the vacuum simulations p is visualized by red isosurfaces on two different contour levels, which enclose 60 and 90% of the probability density (solid and transparent surfaces, respectively).

Display large version of this figure

Figure 8

Transitions of the dye between two subconformations “up” and “down” of conformation A on the nanosecond timescale. The restriction of the mobility of the dye in both conformations is quite similar.

Unfortunately, the occupancy of conformations A and B cannot be directly inferred from the simulation, because only one transition was observed. As a substitute for this lack of reversibility during the simulation time, high-temperature vacuum simulations were carried out to sample the conformational space of the dye more efficiently. From these simulations, a room temperature probability distribution p of the dye positions based on the electrostatic dye-protein interactions was derived, as described in Methods. Figure 7C shows two isosurfaces of p in red enclosing 60 and 90% of the probability density, respectively. Qualitatively, as can be seen from the figure, the dye is in conformation A more often than in B. This is mainly due to the favorable electrostatic interactions between the dye and the protein: In conformation A, the negatively charged dye has contact to two lysines (colored blue in Figure 7C), whereas in conformation B it only has contact to one lysine. In addition, as will be described further below, the calculated fluorescence anisotropy of the dye in conformation A agrees much better with the measured one than it does in conformation B. Thus, we assume that conformation A is the dominant conformation of the dye.

To address the question if and to what extent the dye influences the protein conformation and dynamics, we compared a 5-ns part of the simulation described above to a 5-ns simulation of the same system without the dye. We focused on the change in the flexibility of the protein and in particular of the loop region, where the dye is attached. The flexibility is quantified by the rmsf of the backbone, shown in Fig. 9 for the protein with bound dye (dotted line) and without the dye (solid line). For the calculation of the rmsf the trajectory was fitted onto a reference structure using only the helical residues to focus the analysis on the loop region (thick solid bar). As can be seen, the overall shape of both curves is quite similar. Only the loop residues show a slight decrease of the flexibility, whereas the rest of the protein is hardly affected by the dye. But also for the loop region, the attached dye decreases the rmsf only by maximally 15%. Therefore, the assumption necessarily made for all spectroscopic measurements of this kind, that the dye does not severely influence protein dynamics is, at least for this case, justified.

Display large version of this figure

Figure 9

Root mean square fluctuations (rmsf) of the protein backbone for the protein with the bound dye (dotted line) and without the dye (solid line). The thick solid bar on top denotes the loop residues.

With these results at hand, we can now address the question of which protein region or, more precisely, the motion of which protein region is actually probed by the dye. To this aim, we quantified the correlation between protein and dye motion, as described in Methods and shown color-coded in Fig. 10. Those residues that show the largest correlations (red) are, therefore, those which are “seen” via the dye. Note that the obtained measure of correlation is only a relative one.

Display large version of this figure

Figure 10

One Alexa488 dye molecule attached to the loop of the bR fragment. The protein is colored according to the relative correlation (arbitrary units) of its motion with the motion of the dye. High (low) correlations are shown in red (blue). The bottom panel shows the relative correlation (black line) and the mean square fluctuation (msf; blue line) of the protein backbone.

The highest correlations are seen for residues Gly-33 and Val-34, which are next to the Cys-35 residue to which the dye is covalently bound (Fig. 10, bottom, black curve). These highly correlated residues are in close vicinity to the dye headgroup and also form frequent sterical contacts to the dye. Accordingly, they transmit their high flexibility via noncovalent interactions to this headgroup containing the chromophore of the dye molecule, which is observed in the experiment. Comparison of this correlation with the mean square fluctuation (msf) of the involved residues (Fig. 10, bottom, blue curve), shows that the motion of the dye is mainly governed by the more flexible part of the loop, residues 33 and 34. This agrees with the general expectation that mainly the more flexible residues contribute to the fluorescence depolarization. From this follows that depolarization experiments do not necessarily probe only the residue to which the dye is attached, but rather the flexibility of an entire segment—in our case, the flexible region comprising residues 33,34,35,36,37. Also, the region around residue Ile-45, in the middle part of the left helix (helix B) in Fig. 10, shows a significant correlation, whereas the helix connecting this residue to the loop does not show any correlation. This is remarkable and shows how correlated motions can be transmitted via larger helix parts, e.g., necessary for signal transduction in GPCRs. Here, relatively small correlated motions of the helix induce larger correlated motions in residue 45.

One main goal of this work, as mentioned above, was to identify the kinetic component in the anisotropy decay r(t) of the dye that is governed by and, hence, yields information about, the protein dynamics. To this aim, the reorientational dynamics of the transition dipole moment of the dye was studied. To correlate it with the protein dynamics studied above, the orientational dynamics was analyzed and quantified as a function of dye position, as described in Methods. Figure 11A shows the obtained map, viewed from the bottom side with the protein backbone superimposed; the color encodes the contribution to the depolarization decay at 50ps. This period was chosen as a tradeoff between spatial resolution—the dye should not move too far during this period—and sufficient reorientation.

Display large version of this figure

Figure 11

Analysis of the timescale of the dye-protein correlation. The inset shows the viewing direction used in panels A and B. The ribbons in the foreground depict the protein backbone looking upwards from the bottom of the protein. The backbone is colored according to the correlation of the protein motion with the dye motion, as in Fig. 10. The “cloud” in the background shows all the positions that are visited by the center of mass of the dye headgroup during the simulation. The coloring of this “cloud” denotes the contribution to the depolarization via its orientational mobility. Red, green, and blue indicate high, mid-, and low mobility of the dye, respectively. This coloring has been calculated from the original (A) and from a smoothed (B) trajectory.

As can be seen, is nearly uniform throughout the region covered by the dye. Here, the fast reorientational dynamics of the dye dominate possible—and presumably slower—components that could correlate with the protein motions. These have therefore been filtered in Figure 11B, which was derived similarly as Figure 11A, but using a smoothed trajectory, for which fast fluctuations have been suppressed. The smoothing filter has been chosen to suppress the fastest decay (120ps) within the decay curve calculated from the trajectory as described in Methods; thus, the decay curve calculated from the smoothed trajectory lacks the 120-ps component, but still contains components slower than 300ps.

Indeed, the obtained map shows high orientational mobility (red) near the loop region with a maximum close to residues 33,34, which have been identified above to move strongly correlated with the dye. In contrast, for distant dye positions, its orientational dynamics is less pronounced (blue) and, therefore, contributes little to the observed depolarization. We conclude that there is a slow component of the dye dynamics, which is correlated with the protein motion and which is slower than 300ps.

In the measured curve, shown in Figure 12A, three decay components can be resolved: two faster components (ϕ₁=296ps and ϕ₂=805ps) are assumed to arise from the dye motion relative to the protein (8), where the slower of these two components is suggested to be influenced by the local protein flexibility, as analyzed in the previous section. The third (slowest) component of 5ns originates from the overall tumbling motion of the protein-dye complex. To visualize these multiple decay components, the inset in Figure 12A shows a logarithmic plot of the anisotropy decay and the corresponding fit curve. The measured anisotropy decay curve has been additionally fitted using the cone-in-a-cone model (Eq. (6)); the resulting parameters are shown in Table 1 (Experiment).

Display large version of this figure

Figure 12

Anisotropy decay curves from simulation and experiment, normalized by the initial anisotropy r₀. (A) Measured anisotropy decay (convoluted with the instrument response function). The inset shows a logarithmic plot of this curve together with the residuals of the fit to Eq. (6). (B) Anisotropy of the dye in conformation A (solid line) and in conformation B (short dashes) calculated from the MD trajectory. Also shown (long dashes) is the fit curve of the cone-in-a-cone model to the experimental curve, using the parameters shown in Table 1. (C) Anisotropy in the protein frame for conformation A (solid line) and for conformation B (short dashes), calculated from a trajectory that has been fitted onto a reference structure. The long-dashed line shows the same fit curve as in panel B, except with the global rotational correlation time ϕ_G set to infinity. This curve thus corresponds to the measured anisotropy in the protein frame.

Anisotropy decay curves have been calculated from the simulation, as described in Methods. Separate curves have been calculated for the dye in conformation A and B, shown in Figure 12B as solid and short dashed lines, respectively. As discussed above, for the comparison of the anisotropy with the experiment, the first 16ns of the simulation have been used, during which the dye is in conformation A. Because the overall tumbling motion of the dye-protein system is not of interest here, this motion has been suppressed by least-square fitting all protein structures to a reference structure. The resulting trajectory, thus, describes the dye dynamics in the coordinate frame of the protein. Figure 12C shows the obtained anisotropy decay. Also this decay curve was fitted to the cone-in-a-cone model (Eq. (6)), setting ϕ_G=∞, and thereby also suppressing global rotational diffusion. The decay parameters obtained from this fit (A₁, ϕ₁, A₂, and ϕ₂) are shown in Table 1 (Simulation).

Figure 12C shows that the anisotropy of the dye in conformation B (short dashed line) is significantly higher than in conformation A (solid line), as discussed above. Additionally, the elimination of the overall tumbling allows for a fit using two instead of three exponentials (cf. Eq. (6)), which significantly improves the accuracy of the obtained decay parameters for the local motion of the dye. Because the internal protein and dye dynamics can be assumed to be uncoupled to the overall rotational diffusion, we prefer to compare with experiment the parameters obtained from Figure 12C, i.e., from the more reliable fit.

Considering the systematic acceleration of simulated rotational diffusion, which is, as can be seen, comparable to the acceleration of the free dye analyzed above, the agreement of the fit parameters (Table 1) between experiment and simulation is sufficiently good to assign decay components properly. In particular, the amplitude A₁, which describes the dye-cone angle in the cone-in-a-cone model, matches quite well, indicating that the simulation accurately describes the local wobbling of the dye. The second amplitude A₂ will be analyzed in more detail further below.

The rotational correlation time ϕ_G of the global rotation of the protein was obtained from the experiment, as described above, by fitting the cone-in-a-cone model (Eq. (6)) to the anisotropy decay. From the simulation ϕ_G=3.9ns was obtained via ϕ_G=1/(6D_G) from the rotational diffusion coefficient of the protein D_G=4.3×10⁻⁵ps⁻¹, in good agreement with the experimental value of 5ns.

The limited length of the calculated trajectory (16ns) causes a statistical error for the average (Eq. (2)) in the calculation of the anisotropy decay, and therefore also in the fitted cone-in-a-cone parameters 44. The straightforward approach to estimate this error would be to calculate many similar MD trajectories and to analyze the variance of the decay parameters obtained from each of the trajectories. Unfortunately, the calculation of just one trajectory required already four months of computer time, so this option is not practicable.

Instead, we used an approach, which is principally the same as the one suggested before, but differs in the way the trajectories are obtained. Assuming the dye dynamics to be sufficiently well described by the cone-in-a-cone model, the many required MD simulations are substituted by Brownian dynamics simulations of the transition dipole moment diffusion in the cone-in-a-cone model. From these Brownian dynamics simulations, 230 trajectories were obtained, from which 230 anisotropy decay curves were calculated. By fitting the cone-in-a-cone model, Eq. (6), to these curves, 230 decay parameter sets were obtained. Their variances serve as a measure of the expected statistical errors.

The error of the fast rotational correlation time, Δϕ₁=5ps, is significantly smaller than that of the slow correlation time, Δϕ₂=120ps, because the larger rotational diffusion coefficient yields better sampling of the dye cone. The errors of the amplitudes, ΔA₁=ΔA₂=0.02, imply errors of the semicone angles of Δθ₁=Δθ₂=1°.

In the previous sections, the two fast decay components of the anisotropy decay were tentatively attributed to the local wobbling of the dye on the surface of the protein. It was assumed that the second decay component (A₂, ϕ₂, cf. Table 1) probes the motion of the loop to which the dye is attached.

This assumption is tested here by directly calculating from the simulation the “correct” decay component (A_loop, ϕ_loop) that is due to the loop motion. To this aim, the anisotropy decay curve r_pf(t) calculated above in the coordinate frame of the protein (Figure 12C and Figure 13C, dashed line) is compared to the decay curve r_lf(t) calculated in the loop frame (Fig. 13, solid line). The difference between these two curves is only due to the loop motion.

Display large version of this figure

Figure 13

Calculated anisotropy decay of the dye in conformation A in the protein frame (dashed line) and in the loop frame (solid line). The inset shows several superimposed protein structures from the trajectory, which is fitted onto the loop residues.

The motion of the dye in the protein frame can be described as the motion of the dye in the loop frame superimposed by the motion of the loop. In the framework of the cone-in-a-cone model (see Methods), the corresponding anisotropy r_pf(t) in the protein frame is, thus, given by

(17)

To obtain r_lf(t), all structures from the original MD trajectory have been fitted onto the loop residues 30–42. The inset in Fig. 13 shows several such protein snapshots superimposed. The parameters A_loop=0.77 and ϕ_loop=1370ps are determined from Eq. (17) by fitting. The decay time ϕ_loop of this component matches well to the decay time ϕ₂, but the amplitude A_loop is significantly larger than A₂, i.e., the corresponding semicone angle is smaller. From this we conclude, that the loop motion alone is not sufficient to account for the second decay component (A₂, ϕ₂). This indicates that there must be a further decay component in r_lf(t) that additionally contributes to the depolarization at the timescale of ∼1ns.

Because we observed a flipping of the dye between two orientations in conformation A, as described above (cf. Fig. 8), which occurs roughly on the same timescale (three orientation flip events in 10ns) as the slow correlation times ϕ_loop and ϕ₂, we suppose this process to be the missing additional component in r_lf(t). That means, there are actually two processes, contributing to the slow (second) component in the experimental anisotropy, the flipping of the dye orientation and the loop flexibility, which both occur at the same timescale and therefore cannot be resolved by the experiment alone. The straight use of the cone-in-a-cone model to interpret the measured anisotropy decay would therefore overestimate the absolute amplitude of the loop flexibility.

To test whether the wobbling-in-a-cone model is appropriate to describe the local wobbling of the dye, the orientation distribution of the transition dipole moment was calculated for both dye conformations, A and B. Fig. 14 shows these orientation distributions, calculated as described in Methods. Both distributions are unimodal and resemble Gaussian distributions. The width of the distribution in conformation A, 44°, is broader than that in conformation B, 29°. The shapes of the distributions are well approximated by Gaussian distributions. The semicone angles obtained by a wobbling-in-a-cone analysis of the simulated anisotropy, 40 and 27° for A and B, respectively, agree very well to the obtained widths of the distributions, as has already been shown 45. We conclude that the wobbling-in-a-cone model describes the fast local wobbling of the dye sufficiently accurate, if the cone angle is interpreted as the effective width of the orientational distribution of the transition dipole moments.

Display large version of this figure

Figure 14

Orientation distribution of the transition dipole moment of the dye in both conformations A and B, represented by cones placed at the surface of a sphere. The color and the lengths of the cones denote the frequency that the transition dipole moment adopts a certain direction (red, green, and blue means high, mid-, and low frequency, respectively). The two orientation distributions are centered at the mean position of the dye in conformation A or B, respectively.

Aiming at an atomistic interpretation of time-resolved fluorescence anisotropy experiments, we have carried out molecular dynamics simulations in close resemblance to the experimental situation. Here we have focused on the analysis of fluorescence depolarization measurements of the loop mobility of bacteriorhodopsin (bR). In these experiments, an Alexa488 dye has been attached to the loop within the AB fragment of bR. Our simulations revealed two possible conformations of the dye. It was observed that the mobility of the dye differs significantly for these two conformations. Analysis of the electrostatic interactions suggested one of the dye conformations to be the dominant one that, therefore, is mainly observed in the experiment.

By comparison with a second simulation of the protein without the dye, the influence of the dye on the dynamics of the protein was found to be small, the only significant effect being a small decrease (maximum 15%) of the loop fluctuations due to the bound dye. This finding supports the usual assumption made in the experiments that the dye does not severely affect the protein dynamics.

To study, vice versa, how the dye motion is influenced by the protein, correlations between the dye and the protein motion were analyzed in more detail, using the LMLA algorithm. This calculation revealed those residues that affect the dye motion and that are therefore mainly probed in the experiment. This information is crucial for the atomistic interpretation of the experiment and cannot be inferred from experiment alone.

Overall, the agreement between calculated and measured anisotropy is very good. The calculated anisotropy in the protein frame shows two decay components of 120 and 980ps compared to 410 and 920ps in the experiment. The first decay time is attributed to the rotational diffusion of the dye in the solvent (methanol); thus this fast component is due to the local wobbling of the dye. It has been described previously 43 that, generally, rotational diffusion coefficients of small compounds are quite sensitive to force-field effects, hence the discrepancy for the rotational diffusion. Our assignment is not affected by this possible artifact.

The second component of 980ps has initially been attributed to the protein flexibility. To further test this assumption, the anisotropy was calculated in the coordinate frame of the loop, which allowed to directly assess the influence of the loop flexibility onto the anisotropy. It was found that the decay time of the depolarization induced by the loop (1370ps) is indeed close to the second component. Further evidence for the assignment of the rotational correlation time to the loop dynamics is provided by NMR experiments, which show that the backbone N-H vectors are involved in a 1-ns dynamics 37. However, the loop motion alone can only explain part of the second component. We propose the missing additional contribution to the depolarization on the timescale of ∼1ns to be due to the transition of the dye between the two conformational substates (“up” and “down”), observed in our simulations. The differentiation between these two processes (the loop flexibility and the conformational transition of the dye) provided an interpretation that is not accessible to the experiment. Furthermore, because these two processes in the particular dye-loop construct that we have investigated contribute to the measured depolarization on the same timescale, straightforward application of the conventional cone-in-a-cone model to the anisotropy decay curve, would overestimate the cone angle of the protein cone, i.e., the loop flexibility.