Molecular dynamics simulation of energy migration between tryptophan residues in apoflavodoxin

Abstract

Molecular dynamics (MD) simulations over a 30 ns trajectory have been carried out on apoflavodoxin from Azotobacter vinelandii to compare with the published, experimental time-resolved fluorescence anisotropy results of Förster Resonance Energy Transfer (FRET) between the three tryptophan residues. MD analysis of atomic coordinates yielding both the time course of geometric parameters and the time-correlated second-order Legendre polynomial functions reflects immobilization of tryptophans in the protein matrix. However, one tryptophan residue (Trp167) undergoes flip-flop motion on the nanosecond timescale. The simulated time-resolved fluorescence anisotropy of tryptophan residues in apoflavodoxin implying a model of two unidirectional FRET pathways is in very good agreement with the experimental time-resolved fluorescence anisotropy, although the less efficient FRET pathway cannot be resolved and is hidden in the contribution of a slow protein motion.

---

Introduction

Flavodoxins are small bacterial flavoproteins, which function as low-potential electron-transfer proteins.¹ Apoflavodoxin, which is flavodoxin without cofactor FMN, and flavodoxin have virtually identical structures.² Apoflavodoxin from Azotobacter vinelandii has been used as a prototype to investigate denaturant-induced protein (un)folding.^3–5 The three-dimensional structure of the 179-residue protein is characterized by a five-stranded parallel β-sheet surrounded by α-helices.⁶ This alpha–beta parallel topology is one of the most common protein folds. Fluorescence of apoflavodoxin arises mainly from its three tryptophans (Trp74, Trp128 and Trp167). Trp74 is located in α-helix 3, Trp128 is close to β-strand 5a and Trp167 is in α-helix 5 of the protein. Using picosecond-resolved fluorescence anisotropy of wild-type apoflavodoxin and of mutant proteins lacking one or two tryptophans, it has been demonstrated that photo-excited tryptophan residues of apoflavodoxin exchange energy through a Förster-type of dipolar coupling mechanism (Förster Resonance Energy Transfer, FRET).⁷ Although one might expect reversible or bidirectional excitation energy transfer between the three tryptophan pairs,⁸ this turns out not to be true for apoflavodoxin. Of the tryptophans, Trp128 is located in the most polar environment and consequently its fluorescence spectrum is most red-shifted causing unidirectional radiationless energy migration from Trp74 and Trp167 to Trp128. Energy transfer from Trp167 to Trp128, residues that are 6.8 Å apart, leads to a rapid decay of the experimental anisotropy signal with a 50 ps transfer correlation time.⁷ Energy transfer from Trp74 to Trp128 is very slow and therefore outcompeted by the rapid transfer to Trp 167. Although all three tryptophans of apoflavodoxin are excited, excitation energy rapidly migrates to Trp128, which acts as energy sink. For a schematic overview of the relative orientation of the tryptophan residues, possible pathways of energy transfer and relevant rate constants we refer to Fig. 1. The orientation of the transition dipole moment in the plane of the indole ring (¹L_a) has been obtained from experimental and theoretical studies.^9,10

Fig. 1 Orientation of tryptophan residues in apoflavodoxin. Bold arrows represent the ¹L_a emission transition dipole moment in tryptophan, which direction in the indole ring is taken from.⁷ Grey arrows correspond to the two main, unidirectional energy transfer pathways with experimental reciprocal rate constants (k_ac)⁻¹ = 6.9 ns and (k_cb)⁻¹ = 50 ps.⁷ FRET from Trp74 to Trp128 has negligible efficiency because of the long distance. The rate constants are subscripted with a, b and c denoting Trp74, Trp128 and Trp167, respectively.

The purpose of this study is to compare the experimental time-resolved fluorescence anisotropy with the one obtained with molecular dynamics (MD) simulations over 30 ns using the X-ray structure of A. vinelandii apoflavodoxin as starting structure.⁶ Because of computational limitations two decades ago, MD simulations on another bacterial flavodoxin could only be carried out over a 375 ps trajectory.¹¹ Analysis of nanosecond-length MD trajectories enables the investigation of both the time evolution of FRET parameters and the motional dynamics of the three individual tryptophan residues, which is difficult to access experimentally because of the energy transfer among the tryptophan residues. In addition, the observed time-resolved anisotropy can be compared to the one that is calculated from the MD trajectories. The FRET model that is implemented in the MD simulations involves two unidirectional transfer steps from Trp74 to Trp167 and from Trp167 to Trp128 with experimental rate constants obtained previously (Fig. 1). Although the comparison of MD simulations with experimental time-resolved fluorescence anisotropy experiments has been reported, previous studies were focused on the reorientational dynamics of tryptophan in water^12,13 and of a dye covalently attached to a protein.¹⁴

Methods

Details of molecular dynamics simulations

The X-ray crystallographic structure of the Cys69Ala flavodoxin II from Azotobacter vinelandii (PDB entry1YOB)⁶ was used as the initial structure of the apoflavodoxin. The coordinates of protein atoms corresponding to subunit A were retained, while the FMN co-factor was removed. All missing hydrogen atoms were added using the Leap module of the AMBER10 suit of programs.¹⁵ The electroneutrality of the simulated system was treated by adding 13 Na⁺ counterions. The protein was further immersed in a cubic box of 6494 TIP3P water molecules. The added water molecules were preliminary relaxed and the resulting system was subjected to energy minimization using 1000 steps of steepest descent and 2000 steps or 3000 of conjugate gradient. Both solute and solvent were gradually heated from 10 K to 298 K over 50 ps. Consequently, a 30 ns molecular dynamics (MD) simulation was conducted.

MD calculations were carried out using the SANDER module of the AMBER10 software package with two different force fields (FF, see below).^16,17 The system was set up under the periodic boundary condition in the isobaric–isothermal ensemble (NPT) with a constant pressure of 1 atm and temperature at 298 K. The long-range electrostatic interaction was calculated based on the Particle Mesh Ewald method.¹⁸ All bonds involving hydrogen atoms were constrained according to the SHAKE algorithm.¹⁹ A cutoff distance of 10 Å for non-bonded pair interaction was used with an integration time step of 2 fs. Snapshots were collected every 0.5 ps and only the data taken from the simulation time of 10 ns to 30 ns was collected for analysis.

To improve the statistics of the MD simulations, we performed three independent 30 ns runs. For the force field from 2003 (FF03),¹⁶ the starting MD structures of the first two MD simulations have the same energy conformation but are different from that of the third one. To take into account the effect of a different initial velocity, we used the same energy conformation as the initial structure for the first two simulations, whilst to consider the effect of the starting geometry, a different energy conformation was used in the third MD run. The different random sets of starting velocities for each MD run were then assigned. All geometrical and dynamical parameters reported here are the average parameters from the three MD runs. To investigate the effect of different force fields we also used an older force field (FF99SBildn from 1999).¹⁷ In this case, the three MD simulations have different starting energy conformations. The MD starting velocities of these three systems are also randomly defined, depending on the different velocities assigned by the random number generator seed. Again, the reported geometrical parameters are the average of those obtained in the three simulations. Geometrical and dynamical parameters are presented in the ESI† to allow comparison with the results obtained with FF03.

Theoretical background and computational analysis method

Time-correlation functions. The experimental time-resolved fluorescence anisotropy, A_obs(t), is computed from the parallel and perpendicular fluorescence intensity decays (for details see ref. 7):

The theoretical anisotropy decay, A_theor(t), can be described by the time-correlation function of the reorientation of the emission transition dipole moment:²⁰

where the factor 2/5 accounts for the maximum theoretical anisotropy as a result of photoselection, P₂[·] is the second-order Legendre polynomial, [small mu, Greek, vector]_a(0) reflects the direction of the absorption transition moment at time t = 0 and [small mu, Greek, vector]_e(t) the direction of the emission transition moment at time t. Absorption and emission transition moments may be considered separately, because light absorption occurs in about 10⁻¹⁵ s, which is much faster than the process of internal conversion from the optically excited state to the emitting state (∼10⁻¹² s) before emission starts. The brackets, 〈·〉, denote an ensemble average. The results obtained with experiment and simulation can to a good approximation be compared to each other using the following equation:²¹

〈P₂[[small mu, Greek, vector]_a(0)·[small mu, Greek, vector]_e(t)]〉 = P₂(cos δ)·〈P₂[[small mu, Greek, vector](0)·[small mu, Greek, vector](t)]〉 (3)

where δ is the angle between the absorption and emission transition moments. In eqn (3) the subscripts on the vector μ are omitted, since they represent the direction of the emission transition moment. Therefore in molecular dynamics simulations, we assume to follow the emission transition moment in time. The angle δ can be obtained by eqn (5), where A₀ is the limiting anisotropy at t = 0.

We added eqn (5) to emphasize that in tryptophan residues there is always a mixture of ¹L_a and ¹L_b absorption transition moments with strongly overlapping absorption bands making A₀ significantly smaller than 0.4.

For a multi-tryptophan protein there are three possible contributions to the depolarization of the fluorescence: rotation of the whole protein, restricted motion of the individual tryptophan residues and inter-tryptophan radiationless energy transfer.²² FRET as source of depolarization will be treated in the next section. Let us consider here the two possible types of motion. Since restricted flexibility inside a protein is an order of magnitude faster than protein rotation eqn (2) can be factorized in three terms:

where C_prot(t) and C_int(t) are the correlation functions of protein rotation and restricted motion, respectively. To gain insight in the dynamic behavior of the tryptophan side chains, described by C_prot(t) and C_int(t), the second-order Legendre polynomial functions of the three tryptophan residues were calculated from the MD trajectories. Since we are dealing with one protein per simulation, the ensemble average must be substituted by the time average. The ensemble average at time t is given by:where N is the length of the trajectory. In the following sections the right-hand expression of eqn (7) is abbreviated as P₂(t). In addition, is replaced by A₀. For the fluorescence anisotropy of single- and double-tryptophan apoflavodoxins we have obtained A₀ = 0.24 (ref. 7 and 23).

FRET parameters obtained from MD trajectories

The theoretical anisotropy of each tryptophan residue in apoflavodoxin should be treated independently, because they will be excited simultaneously and separately. Note further that Trp128 is located in a polar environment and acts as an energy sink. Therefore we have mainly unidirectional transfer to this residue, when the other two Trps are excited. Let us first consider the FRET rates between the relevant tryptophan pairs. The tryptophan residues are labeled with a (Trp74), b (Trp128) and c (Trp167), respectively (see Fig. 1). FRET rates strongly depend on geometrical factors like mutual distance (R) and orientation factor (κ²). Subscripted κ² values are the orientation factors between the transition dipole moments of the relevant tryptophan pairs at t = t′ in the trajectory and are given by eqn (8)–(10):

θ^T_ac(t′) is the angle between transition moments of Trp74 and Trp167. θ^R_a(t′) is the angle between the transition moment of Trp74 and the separation vector from the center of Trp74 to the center of Trp167. θ^R_c(t′) is the angle between the transition moment of Trp167 and the same separation vector in opposite direction. The angle θ^T_cb(t′) is the angle between transition moments of Trp167 and Trp128. The other angles have similar meaning as for Trp74 and Trp167, but now for the Trp167–Trp128 couple. Note that in eqn (9) is used, since it is different from the one in eqn (8). For the Trp74–Trp128 couple we have the same definitions for all three angles. Note that in eqn (10) and are used to distinguish them from corresponding angles in eqn (8) and (9).

The FRET rate constants at t = t′ from a to c, from c to b and from a to b are given by eqn (11)–(13), representing the famous Förster equation for the rate of resonance energy transfer:

R^ac₀(t′) is the critical transfer distance between Trp74 and Trp167, R^cb₀(t′) is the one between Trp167 and Trp128 and R^ab₀(t′) is the one between Trp74 and Trp128. All critical transfer distances are a function of the orientation factor in the trajectory at t = t′. R^ac(t′) is the center-to-center distance between Trp74 and Trp167, R^cb(t′) the one between Trp167 and Trp128 and R^ab₀(t′) the one between Trp74 and Trp128, at t = t′. τ_a and τ_c in eqn (11)–(13) are the fluorescence lifetimes of Trp74 and Trp167, respectively, without FRET. τ_a is known from previous work, τ_a = 3.2 ns.^7,23 τ_c is unknown, but can be estimated as follows. The average fluorescence lifetime of the WFW apoflavodoxin mutant has been determined as 2.6 ns.⁷ In first approximation the observed lifetime of 2.6 ns can be considered as a mean value of the lifetimes of both tryptophan residues yielding a fluorescence lifetime τ_c = 2.0 ns.

Analysis of MD anisotropy trajectories including FRET

Let us now consider the anisotropy decays arising from the individual tryptophan residues.

(i) When Trp128 residue is excited, it will not show resonance energy transfer uphill to Trp167 and Trp74. The only depolarization is by protein rotation. The expected anisotropy can be approximated by:

A_bb(t) = A₀P^b₂(t) (14)

where P^b₂(t) is the 2nd order Legendre polynomial for Trp128 that can be evaluated via eqn (7).

(ii) Excitation of Trp167 will lead to immediate energy transfer to the nearby Trp128 with a transfer rate constant k_cb, followed by further depolarization through protein rotation. There will be simply no time for resonance energy transfer to Trp74. The expected anisotropy can be approximated by:

A_cb(t) = 〈[β^T_cb(t′)exp{−k_cb(t′)t} + β^R_cb(t′)]〉_AVP^b₂(t) (15)

〈·〉 is an averaging procedure with respect to t′ over the whole trajectory. The amplitudes β are defined below.

(iii) Excitation of Trp74 will give resonance energy transfer mainly to Trp167 with a transfer rate constant k_ac (unidirectional), after which excitation is immediately transferred to Trp128. The rate constant of transfer is the one for the Trp74–Trp167 couple, but the amplitudes β are for the Trp74–Trp128 couple. The expected anisotropy can be approximated by:

A_ac(t) = 〈[β^T_ab(t′)exp{−k_ac(t′)t} + β^R_ab(t′)]〉_AVP^b₂(t) (16)

with the same averaging procedure as described for eqn (15).

Because of the longer distance, resonance energy transfer from Trp74 to Trp128 is much too slow to compete with transfer from Trp74 to Trp167. However, effectively there is almost direct transfer from Trp74 to Trp128 via the Trp167-detour. The excited state of Trp167 is hardly being populated because of immediate transfer to Trp128 and is, expectedly, not observed in fluorescence because of the very low quantum efficiency. Note further that in all three expressions the correlation function of P^b₂(t) appears because it is the only emitter after FRET.

The amplitudes β in eqn (15) and (16) obey the following relationships:

where all angles are defined in eqn (9) and (10).

In order to compare with the observed anisotropy the calculated anisotropy obtained by MD simulation, A_calc(t), can be expressed by:

A_calc(t) = f₁A_bb(t) + f₂A_cb(t) + f₃A_ac(t) (21)

The fractions f_i (i = 1–3), the sum of which is normalized to unity, are related to the relative absorbances of the tryptophan residues at the excitation wavelength (300 nm) and relative fluorescence quantum efficiencies at the emission wavelength (349 nm). The fractions f_i (i = 1–3) are varied to obtain the minimum value of the quality of fit criterion between the observed and the calculated anisotropies. Fitting was achieved with a Microsoft Excel spreadsheet using the Solver routine. A least absolute values approach (LAV) was used during the analysis.²³

Results

Geometrical fluctuations of tryptophan residues

The atomic coordinates obtained from the MD simulation yield useful FRET information. Since there are three Trp pairs, we can measure the distance fluctuations of the three pairs over the whole simulation period, as well as the fluctuations in the inter-planar angle of each pair. Graphical representations of distance fluctuations are shown in Fig. 2.

Fig. 2 Fluctuating center-to-center distances between tryptophan pairs in apoflavodoxin. Green: Trp74–Trp167; blue: Trp74–Trp128; red: Trp167–Trp128. Average values are collected in Table 1.

Table 1 Geometrical factors of tryptophan pairs in apoflavodoxin obtained from the MD trajectories

Trp pair^a	Distance^b (nm)	Inter-planar angle^c (deg)	κ^2 (—)
a a, b and c denote Trp74, Trp128 and Trp167, respectively.b Distance is center-to-center distance.c Inter-planar angle is between the two indole rings of the pair. The inter-planar angle is not defined for cb and ac.
cb	0.743 ± 0.025	—	0.45 ± 0.14
ac	1.461 ± 0.019	—	0.40 ± 0.11
ab	2.007 ± 0.027	18.6 ± 4.7	0.24 ± 0.12

---

One important conclusion is that the inter-tryptophan distances hardly change over the simulation period of 20 ns.

The inter-planar angles in the three MD simulations are presented in Fig. 3.

Fig. 3 Fluctuating inter-planar (indole) angles between tryptophan pairs in three different MD simulations. Green: Trp74–Trp167 (ac); blue: Trp74–Trp128 (ab); red: Trp167–Trp128 (cb).

The inter-planar angle between Trp74 and Trp128 is constant. Interestingly, Trp167 is flipping to another stable orientation for the inter-planar angles between Trp167 and the other two tryptophan residues. These sudden flip-flops will be discussed separately.

The orientation factors between the transition moments of each pair (see eqn (8)–(10)) exhibit fluctuations between 0 and 1.0 (Fig. 4), although they are apparently fluctuating around a constant average value.

Fig. 4 Fluctuating orientation factors κ² between two different tryptophan pairs in apoflavodoxin. Green: κ²_ac; blue: κ²_ab; red: κ²_cb. The orientation factors are obtained using eqn (8)–(10). Average values are collected in Table 1.

All averaged geometrical factors are collected in Table 1. In conclusion, these simulations show that native apoflavodoxin is a rigid protein, although there is room for a flip of Trp167 on nanosecond timescale.

Fluctuations of rate constants of energy transfer between different pairs

The critical distances have been determined for all three Trp–Trp FRET couples.⁷ The values are R^ac₀ = 1.26 nm, R^cb₀ = 1.42 nm and R^ab₀ = 1.11 nm. The differences are due to different orientation factors κ² for the three pairs as determined from the X-ray structure. We may correct these critical distances with the κ² values listed in Table 1 that are representing MD-averaged orientation factors. The new values are R^ac₀ = 1.24 nm, R^cb₀ = 1.27 nm and R^ab₀ = 1.14 nm. However, we will not use the critical distances like this. The critical distance R₀ is proportional to (κ²n⁻⁴JQ_d)^1/6, where κ² is the orientation factor, n is the refractive index of the medium between the tryptophan residues in the wavelength range of the spectral overlap integral J and Q_d is the fluorescence quantum yield of the donor in the absence of transfer. When we assume that n, J and Q_d are invariant during the simulations, then R₀ only depends on (κ²)^1/6. The critical distance for κ² = 1 can be calculated for each pair amounting to R₀ = 1.45 nm. Now we can multiply this number with (κ²)^1/6 for each time point in the trajectory of the three pairs (see Fig. 4) to obtain the fluctuation in R₀ for each pair. Then, we use these results together with the distance fluctuations (Fig. 2) and the fluorescence lifetimes in the absence of FRET, τ_a (3.2 ns) and τ_c (2.0 ns), to obtain the fluctuating rate constants of resonance energy transfer of each pair with the help of eqn (11)–(13). The fluctuating rate constants in the MD trajectory are presented in Fig. 5. The FRET rate constants averaged over the MD trajectory are: k_cb = 24.7 (ns)⁻¹, k_ac = 0.265 (ns)⁻¹ and k_ab = 0.035 (ns)⁻¹. The distance is the most important factor for the rate constants. The large distance between Trp74 and Trp128 (2.0 nm) nicely illustrates that the FRET rate constant is orders of magnitude smaller and its contribution can be safely ignored.

Fig. 5 Trajectories of rate constants of energy transfer between two different tryptophan pairs in apoflavodoxin plotted on a semilogarithmic scale. FRET rate constants are obtained with eqn (11)–(13). Green: Trp74–Trp167; blue: Trp74–Trp128; red: Trp167–Trp128. The FRET rate constants averaged over the MD trajectory are k_cb = 24.7 ± 5.4 (ns)⁻¹, k_ac = 0.256 ± 0.025 (ns)⁻¹ and k_ab = 0.035 ± 0.005 (ns)⁻¹.

Motional dynamics of tryptophan residues

The MD simulations should also reveal whether the tryptophan residues are rigidly bound to the protein by following the 2nd order Legendre polynomial function, P₂(t), of the tryptophan emission transition moment describing the angular displacement (see eqn (7)). The P₂ correlation functions are presented in Fig. 6.

Fig. 6 Second order Legendre polynomial function, P₂(t), describing the reorientation of the emission transition dipole moment of the three tryptophan side chains in apoflavodoxin calculated with eqn (7). Dotted line: Trp128; solid line: Trp74; dashed line: Trp167.

All three correlation functions decay in the nanosecond time range reflecting motions of large protein segments. In the beginning a small amount of picosecond internal mobility is observed that most likely represents the internal correlation function C_int(t) (eqn (6)) and can be interpreted as ‘rattling in a cage’. The correlation functions indicate rigidly bound tryptophan residues, which information is complementary to that obtained from the geometric parameters. The slow motion of Trp74 is experimentally reflected in the fluorescence anisotropy decay of the single-tryptophan WFF mutant showing that Trp74 rotates together with the whole protein.^7,23 The slowly decaying time-correlation functions are in agreement with NMR relaxation studies, in which tryptophans of other (holo-)flavodoxins are immobilized in the protein matrix.^24,25

Simulation and fitting of time-resolved fluorescence anisotropy

The time-resolved anisotropy calculated from the MD trajectory, A_calc(t), is given by eqn (21) and its three contributing factors are summarized here. A_bb(t) is the contribution to the anisotropy of Trp128 showing no FRET (eqn (14)), A_cb(t) is the contribution to the anisotropy of Trp167 with ultra-rapid FRET to Trp128 (eqn (15); MD-averaged rate constant is 24.7 (ns)⁻¹) and A_ac(t) the contribution to the anisotropy of Trp74 with slow FRET to Trp167, which is followed immediately by rapid FRET to Trp128 yielding amplitudes β_ab instead of β_ac (eqn (16); MD-averaged rate constant is 0.256 (ns)⁻¹). The fractions f_i (i = 1–3), the sum of which is normalized to unity, are weighting factors. The time-dependence of fluctuating amplitudes β are presented in Fig. 7.

Fig. 7 Time-dependence of fluctuating amplitudes β^T and β^R as obtained from the MD trajectory. The relevant pair is indicated above each panel. The solid traces are β^T. The dotted traces are β^R (in panels ac and cb these are the lower curves). The average values and standard deviations are collected in Table 2.

Table 2 Average values of the amplitudes β obtained from the MD trajectory^a

Trp pair	β^T (—)	β^R (—)
a a, b and c denote Trp74, Trp128 and Trp167, respectively. The sum of β^T and β^R is normalized to 0.24.
cb	0.212 ± 0.008	0.028 ± 0.008
ac	0.165 ± 0.024	0.075 ± 0.024
ab	0.119 ± 0.021	0.121 ± 0.018

---

The values averaged over the MD trajectory (eqn (17)–(20)) are collected in Table 2.

We used the average β values listed in Table 2 for the calculation of eqn (15) and (16) (β^T_cb = 0.212; β^R_cb = 0.028; β^T_ab = 0.119; β^R_ab = 0.121) and the averaged transfer rate constants obtained in the previous section. The normalized (to A₀ = 0.24), individual tryptophan contributions to the time-resolved anisotropy as calculated from the MD trajectories are shown in Fig. 8.

Fig. 8 Normalized (to A₀ = 0.24) individual time-dependent anisotropy curves calculated from MD trajectories using eqn (14)–(16). Dashed curve: emission from Trp128 upon excitation of Trp128 (A_bb); dotted curve: emission from Trp128 upon excitation of Trp74 (A_ac); solid curve: emission from Trp128 upon excitation of Trp167 (A_cb).

A_bb possesses the slowest decay because of the absence of FRET, while the decay of A_cb is ultra-rapid because of very efficient FRET. A_ac incorporates the slow transfer step and decays only slightly faster than A_bb. In all contributions the 2nd order Legendre correlation function of Trp128 is used, because that residue is the emitting species after FRET.

The observed and calculated anisotropies are presented in Fig. 9. When the calculated anisotropy is fitted to the observed anisotropy with the fractions f_i as fit parameters (see eqn (21)), the fit is very good with relative fraction f₁ = 0.08 for A_bb(t), relative fraction f₂ = 0.69 for A_cb(t) and relative fraction f₃ = 0.23 for A_ac(t). The values of the fitted fractions will be discussed separately.

Fig. 9 Experimental anisotropy (dots) and calculated anisotropy (solid line) consisting of a linear combination of individual curves in Fig. 8 with fractions f_i (i = 1–3) (eqn (21)) optimized in the fit procedure. Results are: f₁ = 0.08 (A_bb), f₂ = 0.69 (A_cb) and f₃ = 0.23 (A_ac).

Discussion

The MD simulations over a 30 ns trajectory have led to new insight in the molecular properties of the tryptophan residues in apoflavodoxin. The temporal fluctuation of distances and inter-planar angles and the time course of the second-order Legendre polynomial function reflect relative immobilization of two tryptophans (Trp74 and Trp128) in the protein matrix, while there is room for flip-flop motion of Trp167. The structural characteristics of Azotobacter vinelandii apoflavodoxin have been determined using multidimensional NMR spectroscopy.² Apoflavodoxin has a stable, well-ordered core but its flavin binding region is flexible. This dynamic disorder in the flavin binding region is perhaps related to the flip-flop motion. From the data presented in Fig. 3 we can determine the frequency of flip-flop events. For the three 20 ns trajectories of the inter-planar angle between Trp128 and Trp167 we counted on the average 11 events giving a flip-flop frequency of 0.55 (ns)⁻¹. Similarly, for the Trp74–Trp167 pair we counted 13 events yielding a flip-flop frequency of 0.65 (ns)⁻¹. The duration of flip-flop events shows variation and is in the range of 200–800 ps. Although we have only focussed attention to the tryptophan residues, other aromatic residues in the flavin binding region may also exhibit similar dynamic behavior. A good control for the absence of flip-flop motion would be to carry out MD simulations of flavodoxin with the cofactor FMN bound to it.

In the following paragraphs we discuss the significance of the found fractions f₁, f₂, and f₃. The simulated time-resolved fluorescence anisotropy of tryptophan residues in apoflavodoxin invoking two unidirectional FRET pathways (from Trp74 to Trp167 and from Trp167 to Trp128, see Fig. 1) is in very good agreement with the experimental time-resolved fluorescence anisotropy (Fig. 9). The values of optimized fractions are f₁ = 0.08 (A_bb), f₂ = 0.69 (A_cb) and f₃ = 0.23 (A_ac). The energy-transfer efficiency E for the Trp74–Trp167 pair can be easily calculated from:

After substituting MD-averaged distances (R^ac = 1.46 nm; R^ac₀ = 1.24 nm) in eqn (22), the result is E = 0.273 or 27%. This quantity can be determined in practice, as has been shown from previous fluorescence anisotropy decay analysis of wild type and mutant apoflavodoxins.⁷ Since the transfer efficiency amounts to 27%, Trp74 molecules that do not participate in FRET will show depolarization by rotational motion. Let us consider the slow rotational motion of Trp74 in the simulated anisotropy, which must have a similar form as in eqn (14):

A_aa(t) = A₀P^a₂(t) (23)

When A_aa(t), A_bb(t) and A_ac(t) are plotted, we can see that all anisotropy contributions have similar decays in the first four nanoseconds and can not be distinguished from each other in practice (Fig. 10).

Fig. 10 Normalized (to A₀ = 0.24), individual time-dependent anisotropy curves calculated from MD trajectories using eqn (14), (16) and (23). Dashed curve: emission from Trp128 upon excitation of Trp128 (A_bb); dotted curve: emission from Trp128 upon excitation of Trp74 (A_ac); solid curve: emission from Trp74 upon excitation of Trp74 (A_aa).

The noise in the experimental anisotropy decay increases over time (Fig. 9), which can be understood by the decreasing fluorescence intensity (the fluorescence lifetime of wild type apoflavodoxin is 4.0 ns). Thus, the following conclusion is valid to interpret our MD simulations: the less efficient transfer cannot be distinguished from slow rotational motion in the simulated anisotropy decay. This conclusion can be corroborated by two more fitting results constraining either f₁ or f₃ to be equal to zero leaving only two fractions to be fitted. In the first case we obtained f₂ = 0.67 and f₃ = 0.33 and in the second case f₁ = 0.30 and f₂ = 0.70 with similar fit qualities based on the sum of LAV values (20.7 for unconstrained fitting and 20.8 and 20.9 for constrained fitting, respectively).

What is the reason for the appearance of the relatively high fraction f₂(A_cb)? To answer this question it is relevant to discuss some features of the experimental anisotropy decay. The used excitation wavelength was 300 nm, which is at the red edge of the light absorption spectrum. This excitation wavelength was chosen because of the much higher initial anisotropy than with shorter ones. It is very likely that Trp167 possesses the highest absorption cross section at this wavelength. Experimental evidence comes from fluorescence excitation spectra of two available single-tryptophan apoflavodoxin mutants, namely WFF (Trp128 and Trp167 replaced by phenylalanine) and FFW (Trp74 and Trp128 replaced by phenylalanine). A stable FWF mutant apoflavodoxin could not be prepared. FFW apoflavodoxin indeed shows a twice higher absorption cross section than the one of WFF at 300 nm (see Fig. S11†).

Another experimental point to emphasize is the following. The ultra-short transfer correlation time of 50 ps is accurately recovered from analysis of the time-resolved anisotropy.⁷ However, the measured instrumental response function is of the same order of magnitude (∼80 ps full width at half maximum). This leads to distortion of the anisotropy decay and, after deconvolution, to lower initial anisotropy values.²⁶ Indeed, the experimental A₀ = 0.19 is obtained and not the expected A₀ = 0.24 that is observed with single- or double-tryptophan mutants of apoflavodoxin.⁷ The difference of initial anisotropies may be another source of uncertainty in the determination of the fractions. Therefore, to compare MD simulations with anisotropy experiments exhibiting processes of picosecond duration, one should rely on a technique with sub-picosecond time resolution such as fluorescence up-conversion spectroscopy for which correct initial anisotropies have been obtained.²⁷

Finally, the results with the other force field (FF99SBildn, see Supplementary Information) can be summarized as follows. The geometrical and FRET parameters of the three pairs (distance, dihedral angles, orientation factor and transfer rates) are the same (Fig. S1–S5 and Table S1†). The quantities needed for simulation of the time-dependence of the fluorescence anisotropy show some changes (Fig. S6–S10 and Table S2†). The Legendre polynomials of Trp128 and Trp74, in particular, exhibit a significantly faster decay (Fig. S12†), which ultimately results in a poorer fit at longer times (>3 ns) between simulated and observed anisotropies (Fig. S9†).

Conclusion

By performing MD simulations over a similar timescale as fluorescence anisotropy decay measurements, it turns out that the anisotropy model of two unidirectional FRET steps from two tryptophan residues to a third one acting as a sink, can be exactly reproduced from the atomic coordinates of the three tryptophan residues in the MD trajectory. Since the less efficient FRET pathway and protein motion take place on similar timescales, these contributions to the simulated fluorescence anisotropy decay cannot be distinguished.

Acknowledgements

The Supramolecular Chemistry Research Unit, Mahasarakham University, and Computational Chemistry Unit Cell, Chulalongkorn University, are acknowledged for providing computing facilities. The Ratchadaphiseksomphot Endowment Fund of Chulalongkorn University (RES560530217-HR) is also acknowledged. This work was supported in part by the Netherlands Organization for Scientific Research (NWO), Molecule to Cell programme (N.V.V.).

References

1. S. G. Mayhew and M. L. Ludwig, in The Enzymes, ed. P. D. Boyer, Academic Press, New York, 1975, vol. 12B, pp. 57–118.
2. E. Steensma and C. P. van Mierlo, J. Mol. Biol., 1998, 282, 653–666.
3. Y. J. Bollen, I. E. Sanchez and C. P. van Mierlo, Biochemistry, 2004, 43, 10475–10489.
4. Y. J. Bollen, M. B. Kamphuis and C. P. van Mierlo, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 4095–4100.
5. S. M. Nabuurs, A. H. Westphal and C. P. van Mierlo, J. Am. Chem. Soc., 2009, 131, 2739–2746.
6. S. Alagaratnam, G. van Pouderoyen, T. Pijning, B. W. Dijkstra, D. Cavazzini, G. L. Rossi, W. M. Van Dongen, C. P. van Mierlo, W. J. van Berkel and G. W. Canters, Protein Sci., 2005, 14, 2284–2295.
7. N. V. Visser, A. H. Westphal, A. van Hoek, C. P. van Mierlo, A. J. Visser and H. van Amerongen, Biophys. J., 2008, 95, 2462–2469.
8. F. Tanaka and N. Mataga, Photochem. Photobiol., 1979, 29, 1091–1097.
9. B. Albinsson and B. Norden, J. Phys. Chem., 1992, 96, 6204–6212.
10. P. R. Callis, Methods Enzymol., 1997, 278, 113–150.
11. R. Leenders, W. F. van Gunsteren, H. J. C. Berendsen and A. J. W. G. Visser, Biophys. J., 1994, 66, 634–645.
12. X. Daura, R. Suter and W. F. van Gunsteren, J. Chem. Phys., 1999, 110, 3049–3055.
13. P. Mark and L. Nilsson, J. Phys. Chem. B, 2002, 106, 9440–9445.
14. G. F. Schröder, U. Alexiev and H. Grubmüller, Biophys. J., 2005, 89, 3757–3770.
15. D. A. Case, T. A. Darden, T. E. Cheatham, C. L. Simmerling, J. Wang, R. E. Duke, R. Luo, M. Crowley, R. C. Walker, W. Zhang, K. M. Merz, B. Wang, S. Hayik, A. Roitberg, G. Seabra, I. Kolossvary, K. F. Wong, F. Paesani, J. Vanicek, X. Wu, S. R. Brozell, T. Steinbrecher, H. Gohlke, L. Yang, C. Tan, J. Mongan, V. Hornak, G. Cui, D. H. Matthews, M. G. Seetin, C. Sagui, V. Babin and P. A. Kollman, AMBER 10, University of California, San Francisco, 2008.
16. Y. Duan, C. Wu, S. Chowdhury, M. C. Lee, G. Xiong, W. Zhang, R. Yang, P. Cieplak, R. Luo, T. Lee, J. Caldwell, J. Wang and P. Kollman, J. Comput. Chem., 2003, 24, 1999–2012.
17. K. Lindorff-Larsen, S. Piana, K. Palmo, P. Maragakis, J. L. Klepeis, R. O. Dror and D. E. Shaw, Proteins, 2010, 78, 1950–1958.
18. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
19. J. P. Ryckaert, G. Cicotti and H. J. C. Berendsen, J. Comput. Phys., 1997, 23, 327–341.
20. T. Ichye and M. Karplus, Biochemistry, 1983, 22, 2884–2893.
21. P. H. Axelsen, C. Haydock and F. G. Prendergast, Biophys. J., 1988, 54, 249–258.
22. F. Tanaka, J. Chem. Phys., 1998, 109, 1084–1092.
23. S. P. Laptenok, N. V. Visser, R. Engel, A. H. Westphal, A. van Hoek, C. P. van Mierlo, I. H. van Stokkum, H. van Amerongen and A. J. Visser, Biochemistry, 2011, 50, 3441–3450.
24. A. Hrovat, M. Blumel, F. Lohr, S. G. Mayhew and H. Rüterjans, J. Biomol. NMR, 1997, 10, 53–62.
25. P. Zhang, K. T. Dayie and G. Wagner, J. Mol. Biol., 1997, 272, 443–455.
26. J. Papenhuijzen and A. J. W. G. Visser, Biophys. Chem., 1983, 17, 57–65.
27. X. Shen and J. R. Knutson, J. Phys. Chem. B, 2001, 105, 6260–6265.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c4ra03779k

---