INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS FLUORESCENCE DECAYS SIMULATION RESULTS DISCUSSION APPENDIX REFERENCES

Fluorescence intensity and anisotropy decays of proteins from low concentration samples can be accurately measured using modern techniques provided they contain an aromatic amino acid such as tryptophan (Trp), tyrosine, or phenylalanine. However, until recently, the data could not be easily interpreted in part due to their extraordinary sensitivity to the environment. Consequently, fluorescence techniques have contributed much less to structural determinations than molecular modeling, NMR, let alone x-ray crystallography. When it is possible to predict fluorescence properties from a conformation, back-calculations can be iterated, finally yielding a refined insight into the structures and dynamics of proteins in solution.

X-ray, time-resolved fluorescence and molecular dynamics (MD) simulation have recently been successfully combined to study the Human Immunodeficiency Virus-1 (HIV-1), protease/inhibitor complex (Ringhofer et al., 1999). Binding of the inhibitor induced a faster decay of both the experimental and computed protease fluorescence anisotropy decays. The total anisotropy decay of the four Trp system was assumed to be a linear combination of individual fluorescence anisotropy decays that were deduced from the trajectories of each Trp transition moment according to the pioneer work of Ichiye and Karplus (1983). There was no attempt to take into account excitation migration among Trp residues, probably because transfer rates must be evaluated, and these, at their turn, depend on fluorescence lifetimes.

In this paper, both intensity and anisotropy decays are simulated, which implies that most of Trp de-excitation processes are evaluated. This is now possible thanks to recent contributions. First, eight amino-acid side chains may act as quenchers of the Trp fluorescence (Chen and Barkley, 1998). Second, the peptide bond itself was recognized as a quencher of the indole fluorescence (Chen et al., 1996). The electron transfer rate essentially depends on the proximity of the indole ring CE3 atom to the carbon atom of the peptide bond carbonyl (Sillen et al., 2000). This role was delineated by comparison of experimental lifetimes and short MD simulations of a unique Trp-containing protein. In the present work, the method is applied to a protein containing four Trp residues, namely the HIV-1 integrase (IN) catalytic core. At each step of the dynamics, the instantaneous fluorescence lifetime was computed. Moreover, the probability that a chromophore remained excited during a given time was computed separately for each absorption time, then averaged over all possible absorption times. Resonance energy transfer makes the simulation more complex and was introduced in the calculations to predict not only the fluorescence intensity decays but also the fluorescence anisotropy decays from the trajectories.

The integration of a proviral cDNA into host DNA is a critical step in the life cycle of the HIV-1 because it ensures expression and perpetuation of the viral genome (Sakai et al., 1993). This essential reaction is catalyzed by the viral enzyme IN that has been shown to be necessary and sufficient for the integration reaction in vitro (Bushman et al., 1990; Brown, 1990). The integration function is composed of two steps, both involving the nucleophilic attack of a phosphoester bond by the lone pair of a hydroxyl group. In the first step, called processing reaction, IN removes two 3' nucleotides from each strand of the linear viral DNA, resulting in overhanging CA ends. In the second step, called strand transfer reaction, the newly formed 3' OH act as nucleophilic agents and attack phosphoester bonds on the opposite strands of the target DNA (Brown, 1997). Recombinant HIV-1 IN, produced in Escherichia coli, can carry out both reactions in vitro in the presence of divalent ions such as Mg²⁺ and Mn²⁺ (Bushman and Craigie, 1991). In the same environment, it can also carry out the disintegration reaction that is the apparent reversal of the transfer step if presented with a synthetic dumbbell-shaped oligonucleotide (Chow et al., 1992).

Three distinct regions have been identified in the HIV-1 IN (288 residues) (Andrake and Skalka, 1996). The N-terminal domain, residues 1-49, contains a conserved HHCC motif that binds zinc in a 1:1 stoichiometry (Zheng et al., 1996). Zinc binding is believed to stimulate the multimerization process that enhances the activity (Lee et al., 1997; Deprez et al., 2000). The central catalytic core domain, residues 50-212, contains the catalytic site characterized by the essential D,D(35)E motif (D64, D116, and E152). Although all three domains are strictly required for processing and strand-transfer reactions, the core domain by itself can catalyze the disintegration step. The C-terminal domain, residues 212-288, contributes to DNA binding in a nonspecific manner and to the oligomerization that is necessary for the integration process (Brown, 1997). The 3D structure of each of the three domains is known at the atomic level. In the absence of structural data on the entire protein, the knowledge of the core domain structure and dynamics remains essential for inhibitor studies by docking techniques. Several x-ray crystallographic structures for the soluble catalytic core containing a mutation at the F185 position in the absence of divalent ions are now known. In the first available complete structure, a flexible loop 139-153 dangled out of the protein (Bujacz et al., 1996). The next structures were more compact (Maignan et al., 1998; Goldgur et al., 1998) and did not change much when the crystal was soaked with divalent ions. The structure of the catalytic core in the crystal of the bidomain sequence 52-288 is hardly modified (Chen et al., 2000).

The dynamics of the HIV-1 IN catalytic core has been simulated in the absence or presence of Mg²⁺ from a crystal structure lacking the flexible loop (Weber et al., 1998; Lins et al., 1999). The flexible loop was then built by analogy with homologous sequence in Rous sarcoma virus. The dynamics simulations of the entire hydrated catalytic core showed that the flexible loop retains some secondary structure. The binding of a second divalent ion does not decrease the flexibility in the region of residues 140-149 (Lins et al., 2000a). The second metal ion is likely to be carried into the HIV-1 IN active site by a DNA strand (Lins et al., 2000b).

In this work, two MD simulations were performed, one in the absence and the other in the presence of the physiological cation Mg²⁺. The starting structures were taken as they can be found in the literature. Fluorescence properties were then deduced from the MD simulations and compared with experimental data.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS FLUORESCENCE DECAYS SIMULATION RESULTS DISCUSSION APPENDIX REFERENCES

(50-212) Catalytic core domain preparation

pET-15b IN^50-212F185K expression vector encoding amino acids 50-212 of mutant soluble HIV-1 IN was generously donated by R. Craigie (Laboratory of Molecular Biology, NIDDK, NIH, Bethesda, MD). His-tagged IN catalytic core protein was overexpressed in E. coli BL21 (DE3) and purified under native conditions essentially as previously described (Jenkins et al., 1996). Briefly, at OD 0.8, bacterial cultures were induced by 1 mM isopropyl beta-D-thiogalactopyranoside and incubated for three hours at 37°C. The cell pellets were resuspended in ice-cold buffer A (20 mM Tris-HCl pH 8, 0.5 M NaCl, 4 mM -mercaptoethanol, 5 mM Imidazol), treated with lysosyme for one hour on ice and sonicated. After centrifugation (30 min, 8000 × g), supernatant was filtered (0.45 µm) and incubated for at least 2 h with 1 ml of Ni-NTA agarose beads (Amersham Pharmacia Biotech, Umea, Sweden). The beads were washed twice with 10 volumes of buffer A, 10 volumes of buffer A + 50 mM imidazol, and 10 volumes of buffer A + 100 mM imidazol. His-tagged proteins were then eluted with buffer A + 1 M imidazol. The proteins were dialyzed overnight against 20 mM Tris-HCl pH 8, 0.5 M NaCl, 4 mM -mercaptoethanol, and 10% glycerol (v/v). Aliquots were rapidly frozen on dry ice and stored at 80°C.

DNA substrates

Oligonucleotide DHIV ^5'TGCTAGTTCTAGCAGGCCCTTGGGCCGGCGCTTGCGCC_3' was purchased from Eurogentec (Liege, Belgium) and further purified on 18% acrylamid denaturing gel. For disintegration assays, 100 pmol of DHIV oligonucleotide was radiolabeled using T4 polynucleotide kinase (New England BioLabs (UK) Ltd., Hitchin, UK) and 50 µCi of [-³²P]ATP (3000 Ci/mmol). Kinase was heat-inactivated, and unincorporated nucleotides were removed by a passage through Sephadex G-10 column (Clontech Laboratories UK, Ltd., Hampshire, UK). NaCl was added to the final concentration of 0.1 M. Radiolabeled oligonucleotides were heated to 90°C for 2 min, and the DNA was annealed by slow cooling to room temperature.

Enzymatic disintegration assays

Disintegration assays were performed in the presence of 0.5 pmol of DHIV and 2 pmol of HIV purified catalytic core respectively (Leh et al., 2000). Enzymatic reactions were incubated at 37°C for 1 h. Products were separated in 15% denaturing polyacrylamid gel and analyzed using a STORM Molecular Dynamics phosphorimager.

Time-resolved fluorescence measurements

The time-resolved emission anisotropy was obtained by recording the two polarized emission decays Ivv(t) and Ivh(t), using the time-correlated single photon counting technique. The excitation light pulse source was a Ti-sapphire subpicosecond laser (Tsunami, Spectra Physics, Mountain View, CA) associated with a third harmonic generator tuned at 299 nm. The repetition of the laser was set down to 4 MHz. The fluorescence emission was detected through a monochromator (SpectraPro 150, ARC) set at 350 nm ( = 15 nm) and a time-correlated single-photon counting card SPC-430 (Becker-Hickl GmbH, Berlin, Germany) was used for the acquisition of both excitation light pulse and fluorescence emission. The function of the instrumental response of the laser pulse (100 ps) was recorded by detecting the light scattered by a water solution. The time scaling was 11 ps per channel and 4096 channels were used. The two polarized components of the fluorescence decay and the instrumental response profile were alternatively collected during 90 and 30 s, respectively, until the total count of the Ivv component reached 22-26 millions (to insure a single-photon counting condition, the counting rate never exceeded 40 kHz). The correction for the monochromator transmission (G-factor = Ivv/Ivh) was determined from N-acetyl-tryptophanamide polarized decays under the same conditions. The microcuvette (volume 50 µl) was thermostated with a Haake type-F3 circulating bath. The catalytic core concentration was 500 nM in a buffer containing 20 mM Tris-HCl (pH 7.2), 0.1% NP-40 (v/v), 150 mM NaCl, 1 mM DTT supplemented with either 10 mM MgCl₂ or 1 mM EDTA. The anisotropy decay parameters were extracted from both parallel Ivv(t) and perpendicular Ivh(t) polarized fluorescence decay components elicited by vertically polarized excitation. The corresponding analysis was performed by the Quantified Maximum Entropy Method (MEM) (Brochon, 1994; Livesey and Brochon, 1987). This gave the distributions h_f() for the fluorescence intensity decay and h_a(_c) for the anisotropy decay shown on Fig. 3 and 4. The experimental intensity decays shown on Fig. 11 were calculated using

(1)

The experimental anisotropy decays shown on Fig. 12 were calculated using

(2)

where r_e(0) is the apparent experimental anisotropy at time 0.

Molecular dynamics

In all simulations, the catalytic core (residues 50-212) of the HIV-1 IN contained the single mutation F185H. The starting structure for the simulation of the dynamics in the presence of one Mg²⁺ ion was taken from the molecule C of the 1bl3 Brookhaven Protein Data Bank file (Maignan et al., 1998) for which the coordinates of all the residues of the catalytic loop are available. The end residues (210, 211, and 212) were added to the structure, which was then minimized using our quasi-Newtonian minimizer (Le Bret et al., 1991). Because the total charge of the protein with the Mg²⁺ cation is +1, one Cl counterion was added to make the elementary box electrically neutral. The 13 water molecules nearest to the Mg²⁺ cation were kept as they are in the PDB file. The starting structure for the simulation of the dynamics in the absence of Mg²⁺ ion was taken from the 2itg PDB file (Bujacz et al., 1996). Residues 50, 211, and 212 were added. In that case, an extra Na⁺ cation was added in the simulation box. The system was embedded in a 50 × 60 × 50 Å³ box of water. The EDIT module of AMBER was modified so that a box of a given size could be filled with TIP3P water molecules having a density of 1 g/ml. For the simulations, there were 3835 (= 3822 + 13) and 3813 water molecules in the presence and in the absence of Mg²⁺ ion, respectively.

The SANDER module of AMBER (P. Kollman, University of California, San Francisco) was used to simulate the dynamics. The AMBER forcefield (Cornell et al., 1995) was complemented by the parameters of Mg²⁺, Na⁺, and Cl ions (Aqvist, 1990). Long-range Coulombic interactions were calculated using the particle-mesh technique (Darden et al., 1993) for Ewald sums (Ewald, 1921) with a cutoff of 10 Å and a grid size of 50 × 60 × 50 Å³. No correction was applied for the neglected long-range Van der Waals interactions because they are expected to be small. Covalent bonds containing a proton were constrained to their equilibrium length using the SHAKE algorithm (Ryckaert et al., 1977). The elementary integration time step was 2 fs. Before dynamics could be reliably produced, the system was prepared by two heating procedures. In the first procedure, only the water molecules were allowed to move, and the temperature was increased from 25 K to 300 K by 25 K steps during 14.6 ps. At the end of the first procedure, water molecules were cooled to 10 K. In a second procedure, the temperature of the whole system was increased from 10 K to 300 K over 18 ps. The velocities were then equilibrated during 13 ps at 300 K before dynamics production started. Frames were recorded at 0.1 ps interval.

Quantum mechanics

Quantum properties were calculated using the ab initio package Gaussian98 (Frisch et al., 1998). Indole, 3-methyl indole and 2-3 dimethyl indole were first optimized in the ground state at the restricted Hartree Fock level of theory using the basis 6  31 + G(d). The ten first transitions of the optical spectrum were then computed using the CIS keyword (configuration interaction of single excited orbitals). Several n-* transitions were found. The first two -* transitions were assigned to the ¹L_a and ¹L_b transitions.

	FLUORESCENCE DECAYS SIMULATION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS FLUORESCENCE DECAYS SIMULATION RESULTS DISCUSSION APPENDIX REFERENCES

A proper simulation of fluorescence decays should follow step by step the detailed history of the protein from the absorption of a photon till the fluorescence emission. In the complete story (Callis, 1997), some events are easier to simulate than others. Our aim was to build a primitive model containing the elements that could be easily simulated. The simulation was then compared to the experimental data. Below, the simplifying hypotheses used to simulate the fluorescence properties of a protein after its dynamics are described:

1.	Although all the aromatic residues of a protein are fluorescent, only the indole ring of Trp residues can absorb light when excitation wavelength is 299 nm. Throughout this work, the only absorbing or emitting species were the indole rings.
2.	The intensity of the exciting laser beam used here was small enough so that single fluorescence photons were counted. At a given time, at most one Trp residue is assumed to be excited.
3.	When the Trp residue is excited, the electronic distribution changes in about the same time as the time interval used here to integrate Newton law (2 fs). However, not all the conformations were recorded. The time interval separating two recorded successive conformations of the system was much longer (0.1 ps). Therefore, a residue was assumed either excited or lying in the ground state. In our simulation, the change occurred instantaneously. Once the chromophore is excited, the electrons are redistributed, which implies that the Coulombic charges, the bond, and valence angle parameters are modified. In a proper simulation, the force field of the excited chromophore should be modified. This concerns the equilibrium and rigidity values for each bond, each valence angle, each dihedral angle even, and the Coulombic charges. In this work, the chromophore was only virtually excited and kept in the simulation the same force field it had in the ground state. This hypothesis is not as natural as the previous ones and is only justified by the fact that MD simulations are, as they stand, time and disk demanding. We simply could not afford to generate a dynamics simulation where the force field changed according to the migration of the excitation.
4.	The dynamics was simulated in a box that had a certain orientation relative to the laboratory. This would introduce privileged orientations because our simulation time was short relative to overall tumbling. We imagined that the same dynamics was replicated in all orientations relatively to the exciting laser beam.
5.	The absorption spectrum was assumed to be the same for all Trp residues in the IN catalytic core whatever their position relatively to the solvent. This simplifying hypothesis is experimentally valid in the case of barnase (Willaert et al., 1992). Because of the last two hypotheses, each Trp residue had the same probability of absorbing light.
6.	When a Trp residue absorbed light at time  of the simulation, it had a certain probability p to be in the excited 1La state and (1  p) to be in the excited 1Lb state. When the excitation wavelength is 299 nm, p is certainly close to 1 (Valeur and Weber, 1977). However, the possibility that some of the residues may be excited in the 1Lb state is worth studying. Besides, such a tool may be useful if the excitation wavelength is modified in future works. The absorbing transition moments ma() and mb() (in 1La and 1Lb states, respectively) were then computed from the coordinates of the Trp atoms and the 1La or 1Lb transition charges (see Table 1 and Fig. 1) according to (3) If the chromophore is excited to the 1Lb state, it should return instantly to the 1La state.

View larger version (24K): [in this window] [in a new window]   FIGURE 1   Optimized ground state geometry of the 3-methylindole with Trp atoms nomenclature. Trp CA carbon atom takes the place of the in-plane hydrogen of the methyl group in 3 of 3-methylindole. The 1La transition moments ground|er|excited for indole and 3-methylindole are drawn as given by the ab initio program Gaussian98 (see Materials and Methods). A small error is introduced in the direction of the transition moments when, in a first step, the Mulliken transition electron distribution is centered on the nuclei as shown on Table 1 and, in a second step, Eq. 3 is used. The moment (not drawn for clarity) lies 10° from the 3-methylindole moment and 4° from the indole moment. The x and y axes mentioned in Table 1 are drawn.

The fluorescence intensity that can be measured at  + t when a chromophore absorbs light at , is proportional to the population that remains excited at  + t. Therefore the history of an excitation for the conformations between  and  + t needs to be followed and then averaged over . Here, four competing de-excitation pathways were considered.

De-excitation through emission of a photon

De-excitation through proton transfers

De-excitation through emission of an electron

(4)

De-excitation through fluorescence resonance energy transfer

(5)

(6)

(7)

(8)

DA

(9)

28

DA

(10)

DA

DA

Computation of the fluorescence intensity decay

The fluorescence intensity at time t is proportional to the amount of chromophores that are still excited, at time t, after absorption. Therefore, the populations X_i (i = 1, ... , N) of excited residues must be kept track of frame after frame. Here, N, is the number of Trp residues and is either 4 or 8 according to the monomeric or dimeric state of the IN catalytic core. The frames were recorded at time interval  (namely 0.1 ps). The conformation of the protein is assumed to be the same between times '  /2 and ' + /2. In that case, the populations X_i follow a system of N linear equations,

(11)

When i and j differ, K_ij is exactly the resonance energy transfer rate k_ij from the donor j to the acceptor i. The diagonal terms depend on the instantaneous fluorescent lifetimes and the sum of resonance energy transfer rates from i to j,

(12)

Because the matrix K is assumed to be constant during '  /2 and ' + /2, setting (see Appendix),

(13)

we have

(14)

The computations should be done on N vectors X, each of them corresponding to a different absorbing chromophore. It is convenient to introduce an N × N matrix, F, to perform all the computations simultaneously. The population of excited residue i at  + t, t = n when j has been excited at , is the element F_ij of the matrix F(,  + n). Because each of the N Trp residues has the same probability of absorbing, the matrix F(, ) is a scalar matrix with diagonal elements equal to 1/N. F(,  + n) is obtained by successive multiplications,

(15)

Finally the simulated fluorescence intensity I(t) is the sum of all elements of the average over  of the matrix F(,  + t),

(16)

The brackets   call for the averaging over . If the total number of frames is L, the averaging for n takes into account L  n values. Therefore, the decay points are the average of more values and are more accurate for smaller values of t (or of n). The decay (Eq. 16) has to be compared to the experimental one, corrected from the shape of the flash, I_e(t) (Eq. 1).

Simulation of the anisotropy decay

The anisotropy at time t is computed as follows:

(17)

In this formula, the Trp residues i, are supposed to be excited at time  with probability p in their ¹L_a band (unit transition moment µ_a) and with the probability (1  p) in their ¹L_b band (unit transition moment µ_b). All residues j emit from their ¹L_a band at time  + t. As above, F_ij(,  + t) is the probability that j is still excited at time  + t when i was excited at time . The expression between brackets   is averaged over the values of . It is more accurate for smaller values of t. r₀ is the fundamental anisotropy. The calculated anisotropy decay (Eq. 17) has to be compared to the experimental one, r_e(t) (Eq. 2).

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS FLUORESCENCE DECAYS SIMULATION RESULTS DISCUSSION APPENDIX REFERENCES

Enzymatic activity

Recently, we prepared an entire HIV-1 IN that catalyzes 3' processing, strand transfer and disintegration reactions in the presence of the likely physiological cation, Mg²⁺ (Leh et al., 2000). The (50-212) core domain of IN was purified in the same condition (i.e., in absence of detergent) and assayed for its activity on a dumbbell-shaped DNA substrate. The activity of the protein was followed by the appearance of a short 14-mer oligonucleotide cleaved from the 38-mer DNA substrate (Chow et al., 1992). Results of this assay in the presence of either Mn²⁺ or Mg²⁺ at increasing temperature are shown in Fig. 2. The (50-212) truncated protein efficiently performed the disintegration reaction using a dumbbell substrate, although it was active only when Mn²⁺, but not Mg²⁺, was present. Although some disintegration product could be detected after 60 mn at 4°C, the reaction yield showed a constant increase over the temperature range reaching a maximum at 37°C.

View larger version (12K): [in this window] [in a new window]   FIGURE 2   Disintegration activity of the HIV-1 IN catalytic core expressed as the percentage of product per dumbbell-shaped substrate as a function of temperature, in the presence of 10 mM Mn2+ (squares) and 10 mM Mg2+ (triangles).

Experimental fluorescence decays

The fluorescence intensity and anisotropy decays of the IN catalytic core have been measured in the presence of either Mg²⁺ or the chelating agent EDTA at both 5°C and 30°C. The decays can be decomposed into a distribution of exponentials as shown in Fig. 3 (at 5°C) and Fig. 4 (30°C). There was small effect of the divalent cation on both intensity and anisotropy decays. Temperature also had little effect on the intensity decays. In contrast, short correlation times (<3 ns) were drastically more populated at 30°C. The longest correlation time remained similar at both temperatures in the absence of Mg²⁺, and was even slightly higher at 30°C than at 5°C in the presence of Mg²⁺. This was unexpected because the longest correlation time can be interpreted as the rotational correlation time, _c of a rigid sphere of volume V. It depends on the solvent viscosity, , the absolute temperature, T, and the Boltzman constant, k, according to the Perrin (1929) equation,

(18)

For aqueous solutions, the ratio /T decreases by a factor 2.065 from 5°C to 30°C. The experimental values at 5°C and 30°C of the long rotational correlation time suggested that the volume has doubled. The catalytic core would be monomeric at low temperatures and dimeric at 30°C (see discussion) at a temperature at which it is enzymatically active.

View larger version (22K): [in this window] [in a new window]   FIGURE 3   Experimental determination of the HIV-1 IN catalytic core lifetimes and correlation times at 5°C. (A) and (C) Lifetimes distributions. (B) and (D) Correlation times distributions. The distributions were recovered by the Quantified Maximum Entropy Method as indicated in Materials and Methods: (A) and (B), in the presence of 1 mM EDTA, (C) and (D) in the presence of 10 mM Mg2+.

View larger version (25K): [in this window] [in a new window]   FIGURE 4   Experimental determination of the HIV-1 IN catalytic core lifetimes and correlation times at 30°C. (A) and (C) Lifetimes distributions. (B) and (D) Correlation times distributions. The distributions were recovered by the Quantified Maximum Entropy Method as indicated in Materials and Methods: (A) and (B) in the presence of 1 mM EDTA, (C) and (D) in the presence of 10 mM Mg2+.

Dynamics of the flexible loop

For each frame of both dynamics (with or without Mg²⁺), the dihedral angles _i = C_i1-N_i-CA_i-C_i and _i = N_i-CA_i-C_i-N_i+1 of the protein were compared with those of the starting structures _i0 and _i0. The    root-mean-squared deviation (RMSD) was calculated at each time of the simulation using

(19)

  RMSD from residue n₁ = 51 to residue n₂ = 211 are plotted in Fig. 5. After a half-nanosecond transient evolution, both the   RMSD reached a plateau showing that the dynamics are stable. If the same analysis was done on the loop (from n₁ = 139 to n₂ = 153), the   RMSD was slightly larger (Fig. 6). Moreover, in the absence of Mg²⁺, a second transition was observed around 3.4 ns (Fig. 6 A) where the   RMSD increased from 60 to 80° and seemed to drop down toward the end of the simulation. In the presence of cation,   RMSD was smaller and retained a value around 50° (Fig. 6 B).

View larger version (16K): [in this window] [in a new window]   FIGURE 5     RMSD for all residues of the HIV-1 IN catalytic core during the MD. (A), In the absence of divalent cation. (B), In the presence of Mg2+.

View larger version (20K): [in this window] [in a new window]   FIGURE 6     RMSD for residues 139-153 of the HIV-1 IN catalytic core during the MD. (A) In the absence of divalent cation. (B) In the presence of Mg2+.

Whether a residue is part of a helix can be determined for each frame of the simulation using the DSSP program (Kabsch and Sander, 1983). Most residues remained part of the same structure throughout the simulation. The flexible loop (139-153) had unique behavior in this respect. The end of the flexible loop can take on the structure of the adjacent 4 helix, which, at its largest extension, comprises residues 147 to 168 in crystal structures. In the absence of cation, no structure was observed in our simulation before 3.4 ns (Fig. 7 A). After that time, the 4 helix extended transitorily down to residue 151 (Fig. 7 A). In the presence of cation, the 4 helix extended permanently till V150, and, as a transient feature, till P145 (Fig. 7 B). The rest of the loop (144-139) was always disordered.

View larger version (16K): [in this window] [in a new window]   FIGURE 7   Time-dependent secondary structure analyses for the helicity of residues 139-153. A dot means that the residue belongs to an  helix as determined by the DSSP program. (A) In the absence of divalent cation. (B) In the presence of Mg2+.

Because it is supposed to play a key role in the catalytic activity of the enzyme, the conformation of the D,D(35)E motif has been thoroughly studied. In the first complete crystallographic structure, the oxygen atoms of E152 pointed away from those of D64 and D116 (Bujacz et al., 1996), in contrast with the conformation they had in another cognate system, the avian sarcoma virus catalytic domain (Bujacz et al., 1995). A synthetic way to monitor the relative oxygen conformation is to compute the interaction energy of E152 distal oxygen atoms and D64 and D116. This interaction is the sum of the electrostatic and Van der Waals contributions to the energy. It increases when the oxygen atoms are closer. In the presence of Mg²⁺, the oxygen atoms were in the expected conformation from the very beginning (Fig. 8 B). In contrast, in the absence of Mg²⁺, the E152 distal oxygen atoms pointed away. However, after 3.4 ns, the oxygen atoms became spontaneously closer (Fig. 8 A). This conformational change occurred simultaneously with the   RMSD transition as shown in Fig. 6 A.

View larger version (22K): [in this window] [in a new window]   FIGURE 8   Interaction energy between the E152 distal oxygen atoms and those of D64 and D116 during the MD. (A) In the absence of divalent cation. (B) In the presence of Mg2+.

Dynamics of tryptophan residues

Because of our interest for fluorescence simulation, some important parameters concerning the Trp residues have been computed. To estimate the number of water molecules close to the Trp residues, the following function was computed,

(20)

where the index i runs over the indole atoms and j runs over the water molecule oxygen atoms. It becomes larger when there are more water molecules and when the latter are closer to the Trp residue. Such a function reflects the induced dipole contribution to the Van der Waals interactions between the Trp residue and the water molecules. It was preferred over counting the water molecules, because counting the molecules depends on an arbitrarily chosen cut-off value. From inspection of the conformations using visualization systems, a value of f about 0.07 corresponds to a Trp residue completely exposed to water. The average value of f and its fluctuation over time are reported for each residue in Table 3. Therefore, during both dynamics, in presence and in absence of Mg²⁺, W131 and W132 were completely exposed to water, W108 was half exposed, and W61 was buried. Therefore, when simulating energy-transfer rates, the fully solvent exposed chromophores, W131 and W132, were treated as poor donors, whereas the less exposed chromophores, W61 and W108, were treated as good donors. Moreover, Table 3 shows in which secondary structure the Trp residues were involved. All four Trp belonged to the same secondary structure during the MD. Figure 9 shows the evolution of the dihedral angles _A = C-CA-CB-CG and _B = CA-CB-CG-CD1 (see Fig. 1 for atom nomenclature) for the Trp side chains in the presence and absence of cation. Both dihedrals of W61 and W108 kept the same values throughout the simulations in contrast with those of W131 and W132.

View larger version (39K): [in this window] [in a new window]   FIGURE 9   Evolution of the side chain dihedral angles A = C-CA-CB-CG (left column) and B = CA-CB-CG-CD1 (right column) for the four Trp residues during the MD. Black dots, in the absence of divalent cation. Gray dots, in the presence of Mg2+. Because gray and black traces largely overlap, black dots are often hidden by gray dots.

Dipolar versus multipolar approximation

The evolution of the instantaneous values of the parameter G is shown on Fig. 10 for each Trp-Trp pair along the MD trajectories in the absence of Mg²⁺. The multipolar approximation (above the diagonal) can be easily compared to the dipolar approximation (under the diagonal). Results were similar in the presence of Mg²⁺ (data not shown). Although both approximations gave similar general trends in most conformations, in rare cases they could be completely different (see W61-W108 around 1 ns). This shows that the dipolar approximation is rather crude. Figure 10 also shows another interesting feature: G changed very rapidly. It could jump from a negative to a positive value within 0.1 ps (see W131-W132 around 4 ns). This implied that it passed through zero in the meantime. The amplitude of the jumps was not constant during the simulation. Most pairs endured several regimes (Fig. 10).

View larger version (38K): [in this window] [in a new window]   FIGURE 10   Time evolution (time unit:ns) of the dimensionless geometric contribution (×104) to the fluorescence resonance energy transfer rates for each Trp-Trp pair displayed in matrix outlay. Above the diagonal, multipolar approximation (G, Eq. 6). Below the diagonal, dipolar approximation (G, Eq. 7).

Fitting the intensity and anisotropy decays

The simulation assumes that the resonance transfer rates remain constant between two recorded conformations. This assumption is not valid as the geometric contribution G jumps from favorable to unfavorable orientations within two recorded frames. Then it becomes necessary to check that the results do not depend on what happens between two successive recorded conformations. Therefore, the dynamics were filled with 10 or 100 conformations chosen randomly among the closest 100 frames. The simulated fluorescence intensity and anisotropy decays were identical to the fourth significant digit (data not shown). Therefore, the rates can be assumed constant (although in actuality they are not).

According to our hypotheses, the fluorescence decay depends on few parameters: _rad, , k₀, and S_DA (see Eqs. 4 and 10). The values for these parameters may first be taken from the literature. The radiative lifetime may easily be estimated for Trp derivatives that have a mono exponential decay. Under these conditions, the radiative lifetime is the ratio of the measured lifetime to the quantum yield. For N-acetyltryptophanamide, _rad is 21.4 ns, and, for N-acetyltryptophan, it is 22.9 ns (Szabo and Rayner, 1980). More recently, an empirical relationship between the radiative lifetime and the wavelength of maximum emission intensity _max was proposed (Sillen et al., 2000): _rad varies from 18.2 ns to 23.5 ns when _max varies from 320 to 345 nm. In the absence of experimental data, _rad was first uniformly set equal to 20 ns for all four chromophores. Sillen et al. (2000) propose for , 1.9 Å¹, and for k₀, the electron transfer rate from the surface of Trp, 25 ns¹. Studies on barnase show that the values of S_DA depend on the environment of the donor (Willaert et al., 1992). Because they are surrounded by water, W131 and W132 are bad donors and their S_DA was set to 133,000 Å⁶/ns. W61 and W108 are considered as good donors and their S_DA was set to 641,000 Å⁶/ns. This set of values gave a steeper fluorescence intensity decay both in the absence (Fig. 11 A, curve 3) and in the presence of divalent cation (Fig. 11 B, curve 3) when compared with the experimental curves obtained at the two different temperatures 5°C and 30°C (curves 1 and 2).

View larger version (16K): [in this window] [in a new window]   FIGURE 11   Comparison of experimental and simulated fluorescence intensity decays. (A) in the absence of Mg2+. (B) In the presence of Mg2+. Curves 1 and 2, experimental decays Ie(t) from Eq. 1 at 5°C and at 30°C, respectively. Curve 3, simulated decay with rad = 20 ns, k0 = 25 ns1,  = 1.9 Å1, favorable SDA = 641,000 Å6/ns and unfavorable SDA = 133,000 Å6/ns. Curve 4, best fit of the experimental intensity decay at 5°C using rad = 20 ns,  = 1.9 Å1, favorable SDA = 6410 Å6/ns and unfavorable SDA = 1330 Å6/ns. (A) curve 4, k0 = 3.7 ns1; (B) curve 4, k0 = 4.6 ns1. (A) and (B) curve 5, same parameters as curve 4, except  = 1.7 Å1. (A) and (B) curve 6, same parameters as curve 4, except  = 2.1 Å1.

To fit the experimental 6-ns fluorescence decay in the absence of cation at 5°C, the electron transfer rate constant k₀ must be set ~3.7 ns¹ and the values of the parameters S_DA divided by 100 (Fig. 11 A, curve 4). In the presence of Mg²⁺, the best fit to the 4-ns decay was obtained when k₀ was set ~4.6 ns¹, and again the values of the parameters S_DA should be divided by 100 (Fig. 11 B, curve 4). At 30°C, a good fit was obtained with k₀ set ~3.1 and 5.2 ns¹ in the absence and in the presence of Mg²⁺, respectively (data not shown). A change of _rad from either 20-22 or 25 ns had no influence on the results (data not shown). Figure 11 also shows simulated decays with  set to 1.7 and 2.1 Å¹ (curves 5 and 6). Clearly, in all cases, the best fit values for , ~1.9 Ź, agreed with Sillen et al. (2000). Apparently, the values of _rad and , to a lesser degree, are not so critical as that k₀. Typically, the best-fitting simulated decays are too slow in the first 200 ps, then too rapid untill about 2 ns, then too slow later on. The optimal value of k₀ depends on the size of the window. If the data are analyzed within 200-ps windows, the electron transfer rate k₀ must be increased to values ~27 ns¹ in closer agreement with the value reported by Sillen et al. (2000).

Once the intensity decay is fitted, the only parameters that can be modified for the anisotropy decay simulation are the population p of ¹L_a excited states and the fundamental anisotropy r₀. Now, for the excitation wavelength of 299 nm, the value of p is close to 1 (Valeur and Weber, 1977). The simulated anisotropy has a very steep drop between 0 and 0.1 ps, in a time that is much less than the time response of the apparatus. Therefore, the simulated point at t = 0 was skipped. The ¹L_a contribution of the simulated anisotropy at 0.1 ps was then scaled to the experimental point to give the curves shown on Fig. 12. The calculated anisotropy decay (Fig. 12, curve 4) obtained using optimized parameters as defined in Fig. 11, was compared to the experimental decays. Because the simulated decays are computed as time-correlation functions from the MD simulations, they are significant only at their beginning, at low values of t. For the 5°C experimental decay in the absence of cation (Fig. 12 A, curve 1), the fit was good. In the presence of cation, the simulated anisotropy was fitted in the same way. Only the fit during the first half-nanosecond was very satisfying (Fig. 12 B; compare curves 1 and 4). The simulated anisotropy decay could not be fitted to the 30°C experimental anisotropy decay with similar parameters (Fig. 12, curve 2). The fit could be obtained if the spectral parameters were drastically increased: S_DA should be set equal to 13,300 Å⁶/ns for a poor donor and to 64,100 Å⁶/ns for a good donor (curve not drawn on Fig. 12 for clarity).

View larger version (12K): [in this window] [in a new window]   FIGURE 12   Comparison of experimental and simulated fluorescence anisotropy decays. (A) In the absence of Mg2+. (B) In the presence of Mg2+. Curves 1 and 2, experimental decays re(t) from Eq. 2 at 5°C and 30°C, respectively. The experimental r(0) were 0.252 at 5°C and 0.223 at 30°C in the absence of Mg2+. The experimental r(0) were 0.237 at 5°C and 0.221 at 30°C in the presence of Mg2+. Curve 3, same as curve 2 with the longest anisotropy correlation times divided by 2. Curve 4, simulated anisotropy decay using rad = 20 ns,  = 1.9 Å1, favorable SDA = 6410 Å6/ns and unfavorable SDA = 1330 Å6/ns. (A) curve 4, k0 = 3.7 ns1 and (B) curve 4, k0 = 4.6 ns1.

As shown before, the catalytic core is probably dimeric at 30°C. It was then interesting to use the dynamics on the monomer unit of the catalytic core to build, at each frame, a chimeric dimer as it is in four crystal structures (Bujacz et al., 1996