This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.105.080051v1 91/7/2647    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Beier, C. Articles by Steinhoff, H.-J. Search for Related Content PubMed PubMed Citation Articles by Beier, C. Articles by Steinhoff, H.-J.

A Structure-Based Simulation Approach for Electron Paramagnetic Resonance Spectra Using Molecular and Stochastic Dynamics Simulations

Fachbereich Physik, Universität Osnabrück, Osnabrück, Germany

Correspondence: Address reprint requests to Heinz-Jürgen Steinhoff, Tel.: 49-541-969-2675; E-mail: hsteinho{at}uos.de.

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSION ACKNOWLEDGEMENTS REFERENCES

 
Electron paramagnetic resonance (EPR) spectroscopy using site-directed spin-labeling is an appropriate technique to analyze the structure and dynamics of flexible protein regions as well as protein-protein interactions under native conditions. The analysis of a set of protein mutants with consecutive spin-label positions leads to the identification of secondary and tertiary structure elements. In the first place, continuous-wave EPR spectra reflect the motional freedom of the spin-label specifically linked to a desired site within the protein. EPR spectra calculations based on molecular dynamics (MD) and stochastic dynamics simulations facilitate verification or refinement of predicted computer-aided models of local protein conformations. The presented spectra simulation algorithm implies a specialized in vacuo MD simulation at 600 K with additional restrictions to sample the entire accessible space of the bound spin-label without large temporal effort. It is shown that the distribution of spin-label orientations obtained from such MD simulations at 600 K agrees well with the extrapolated motion behavior during a long timescale MD at 300 K with explicit water. The following potential-dependent stochastic dynamics simulation combines the MD data about the site-specific orientation probabilities of the spin-label with a realistic rotational diffusion coefficient yielding a set of trajectories, each more than 700 ns long, essential to calculate the EPR spectrum. Analyses of a structural model of the loop between helices E and F of bacteriorhodopsin are illustrated to demonstrate the applicability and potentials of the reported simulation approach. Furthermore, effects on the motional freedom of bound spin-labels induced by solubilization of bacteriorhodopsin with Triton X-100 are examined.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSION ACKNOWLEDGEMENTS REFERENCES

 
Proteins provide a large diversity of functions due to their complex intrinsic structures. The structural hierarchy of these macromolecules comprises distinct internal dynamics that span a wide range of different timescales: Atomic fluctuations, leading to bond-stretching and angle-bending motions, have periods of several femtoseconds, whereas collective motions of associated atoms show correlation times in the picosecond time range. The correlation times of side-chain motions as well as local backbone fluctuations can reach up to a few nanoseconds, whereas rotational correlation times of entire monomeric proteins are in the order of 10^–8–10^–6 s. Finally, conformational switching and associated movements of protein domains generally occur in the nano- to millisecond time range (1–3).

Electron paramagnetic resonance (EPR) spectroscopy is a sensitive tool to detect molecular dynamics in almost all of these time ranges. In the special case of continuous-wave X-band EPR spectroscopy, as utilized here, reorientational correlation times of nitroxide radicals between 100 ps and 200 ns are detectable. In particular, EPR in combination with site-directed spin-labeling (4) is appropriate to study the dynamics of proteins under physiological conditions providing important insights into their functional mechanisms. Spin-label side chains are introduced at selected sites via cysteine substitution mutagenesis followed by modification of the thiol group with a specific paramagnetic nitroxide reagent. To assure unique labeling, accessible native cysteines have to be replaced by similar substitutes like serines. A widely used nitroxide reporter group is (1-oxy-2,2,5,5-tetramethylpyrrolinyl-3-methyl)-methanethiosulfonate (MTS) (5), which reacts with cysteines under formation of a disulfide bond. The bulkiness of the generated side chain, in the following designated as "R1" (Fig. 2), resembles that of a large native amino acid. Nevertheless, both the correct folding and the functional properties of the spin-labeled protein have to be verified experimentally.

View larger version (29K): [in this window] [in a new window]   FIGURE 2  Polar histogram plots of the dihedral angle distributions are assigned to all bonds of the cysteine-bound MTS spin-label (R1) essential for the reorientational motion of R1. The shaded stick-representation shows the structure of R1 (hydrogens neglected). The unpaired electron is located in the -orbital of the nitroxide bond (NX-OX). To estimate the internal flexibility of a bound spin-label during MD simulations under different conditions, important rotameric states of R1 attached to the center of a pentadeca glycine helix (model system) are visualized. The angular scales of all histogram plots are similarly oriented, and a dihedral angle of zero defines the cis (syn-periplanar) configuration. The shaded lines indicate the results of an MD simulation of the model system with explicit water performed at 300 K (50 ns). These plots are compared with the data of two in-vacuo simulations performed at 600 K (black lines, 6 ns) and 300 K (light-shaded solid graphs, 40 ns), respectively.

In principle, a molecular dynamics simulation of a bound spin-label within the specific protein environment provides all required dynamic information to calculate an EPR spectrum. To perform such a simulation, a virtual spin-label has to be introduced at the desired site in a model of the protein structure. However, a direct spectrum calculation, using the method first published by Robinson et al. (8), requires a large set of trajectories, each with a length of several hundred nanoseconds, due to the long transversal relaxation time of the nitroxide. For proteins, the calculation of MD trajectories with the required length is beyond the limits of acceptable technical and temporal expense. A first-generation EPR spectra simulation method was developed by Steinhoff and Hubbell (9), which avoids unrealistic computational effort, nevertheless taking into account the complex motion of a spin-label side chain in a realistic protein environment. This method consists of three steps:

1. An MD simulation to determine an effective potential of the global reorientational motion of the nitroxide.
   
2. Several single particle stochastic dynamics (SD) simulations with consideration of this potential.
   
3. The calculation of the EPR spectrum on the basis of the SD trajectories.
   

The MD simulation is performed at high temperatures (600 K) to focus on the amplitude and anisotropy of the spin-label motion. The required length of the trajectory of atomic positions amounts to a few nanoseconds only, during which the complete accessible conformational space of the nitroxide side chain is covered. Since the essential information is the reorientational motion of the nitroxide, expressed in the Euler angle space, an effective potential energy function is derived from the Euler angle trajectory (,ß,)(t). The final set of reorientation trajectories, each >700-ns long, is calculated by an essentially faster SD simulation algorithm. In the present approach, this is achieved by solving the Langevin equation for reorientational motion of a single particle restricted by the potential energy topology. One-hundred-sixty Euler trajectories, obtained from the single particle simulation, are used to calculate the Larmor frequencies and magnetization trajectories. The EPR absorption spectrum is finally given by the real part of the Fourier transform of the averaged magnetization trajectory.

This article presents an improved simulation pathway with special consideration of MD simulations to justify the sampling at 600 K. The simulation approach is further applied to study the dynamics of nitroxide side chains bound to several positions of the membrane protein bacteriorhodopsin (BR) and to shed light on the relation between the nitroxide dynamics and the binding site structure. BR is a light-driven proton pump located in the cytoplasmic membrane of, e.g., Halobacterium salinarium (10). It consists of seven transmembrane helices enclosing the chromophore retinal (Fig. 1 a). In the native purple membrane, BR molecules are present as trimers, which further aggregate to form a two-dimensional lattice (11). Since the invention of crystallization techniques to promote three-dimensional BR/lipid crystals, the structure of BR has been determined by different groups at adequate resolution. Besides analyses of the MD simulation approach, the presented method for EPR spectra simulation is verified here by comparison of calculated spectra, based on the E-F loop conformation consistently presented in different crystal structures of BR, to those spectra experimentally obtained under physiological and solubilized conditions. For this purpose, the EPR spectra of 18 single mutants with a nitroxide side chain at positions 157–171 (Fig. 1 a) were studied. The EPR experiments, performed at 9.5 GHz (X-band), and the interpretation of the experimental data of these mutants, were described elsewhere (12).

View larger version (70K): [in this window] [in a new window]   FIGURE 1  Structure of the bacteriorhodopsin monomer (1BRR, Chain B) as published by Essen et al. (13). (a) All investigated positions (157–171) are located in the cytoplasmic part of helices E and F or in the connecting loop. Retinal, the linear chromophoric group, is visualized in gray. (b) Superpositions of six different models of bacteriorhodopsin obtained by x-ray and EM techniques (black: 1BRR (13); blue: 1C3W (32); red: 1CWQ (31); gray: 1QHJ (35); yellow: 2AT9 (54); and green: 1FBB (33)).

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSION ACKNOWLEDGEMENTS REFERENCES

View larger version (54K): [in this window] [in a new window]   FIGURE 9  Superposition of the wild-type structure (1BRR, shaded ribbon) (13) and the triple-mutant structure (1FBK, black/dark-shaded ribbon) (33) of BR representing the ground state and a structural model of BR with an outwardly tilted helix F, respectively. Results of the motional freedom of R1 bound to the positions 167 (a), 170 (b), and 171 (c) are presented. The C atoms of all three positions in both structural models are shown in ball-representation. The outward tilt of helix F in 1FBK facilitates an increased motional freedom of R1 at position 167, 170, and the inwardly oriented R1 side chain at position 171. (a,b) The black mesh contours (smoothed) represent almost the entire accessible space (>90%) of the nitroxyl nitrogen atom (NX) during the in-vacuo MD simulations (600 K, 8 ns) of R1 at positions 167 and 170, whereas the solid contours depict states of higher residence probability. The simulated spectra (shaded line) are compared with the measured spectra of the solubilized spin-labeled BR (solid line) on the right-hand side. (c) The three contours in the left figure indicate almost the entire sampled space of NX during MD simulations (600 K, in vacuo, 8 ns) with different starting orientations of R1. The shaded contour refers to the MD simulation of the inwardly oriented spin-label bound to the wild-type structure of BR (see Fig. 5), whereas the solid contour marks the motional freedom of R1 bound to the triple-mutant structure (inwardly oriented), and the solid dotted contour is the result of the outwardly oriented R1 bound to the triple-mutant structure. Due to a high steric barrier between the in- and outward orientations of position 171, separate simulations for both cases were necessary. Simulated EPR spectra obtained from both MD simulations of the triple-mutant structure (inwardly oriented, shaded dotted; outwardly oriented, shaded line) are compared with the measured spectrum (middle position, solid line) of R1 bound to position 171 of the solubilized BR.

1. Simulations with a single neutral-capped R1 amino acid within 755 molecules of water.
   
2. Simulations of a pentadeca glycine -helix with a R1 substitution at position 10 surrounded by 948 water molecules.
   
3. In-vacuo simulations of protein subsets at different temperatures.
   

The used MD simulation conditions of 600 K, in vacuo, and additional restrictions are an outcome of detailed investigations presented in Results and Discussion. The final system temperature was achieved after two heating procedures. A first MD simulation started with atomic velocities assigned randomly from a Maxwell distribution at 400 K. This high initial temperature required the application of position restraints on all atoms to prevent structural bursts. Despite these restrictions, collective vibrational modes were able to establish during a run time of 6 ps. During the following simulation, the temperature was increased from 400 K to 600 K by weak coupling to a 600 K heat reservoir. This final heating step was performed with position restraints assigned only to the backbone atoms N, C, O (see Fig. 2) of all amino acids. Although the desired temperature was generally achieved within the first 15 ps, the applied heating time of 30 ps ensured a sufficient equilibration determined by converging total energy values.

For the main simulations of at least 6-ns length, time steps of 2 fs were used while all bond lengths were kept constant using the SHAKE algorithm (20). Due to the recent gain in computational power, simulation periods of 10–20 ns are recommended (see also Fig. 3 b). To reach rare conformational states and to assure statistical sufficient sampling within 6 ns, all force constants and Lennard-Jones coefficients were kept unchanged. All heavy atoms of the polypeptide main chain except the -carbons (backbone N, C, O) of helices or ß-sheets were restricted by harmonic position restraints. The restraining force constant was chosen to allow a mean deviation of 0.01 nm from the reference position at 600 K. For soluble proteins, weaker position restraints are suggested. If a crystal structure obtained at room temperature is available, realistic force constants can be estimated from the Debye-Waller factors of the backbone atoms. MD simulations (GROMACS) at 600 K with different force constants were applied to analyze the mean-square fluctuations of the backbone atoms and to relate these values to experimental B-factors considering the equation B = 8²x² [Ų]. Applied position restraint force constants of 500 and 100 [kJ mol^–1 nm^–2] yield mean-square fluctuations that correlate with B-factors of 27 and 74 [Ų], respectively. The position restraints are not only necessary to keep the peptide fragments of the protein subsets at their physiological positions but also to inhibit rotations of the whole subset or helical components. The samples utilized for EPR investigation in this article consist of native purple membrane sheets including a high concentration of BR trimers. Hence, rotation of entire BR molecules as well as significant tumbling of complete helices is negligible below 1 µs. The probable spectral contribution of the overall rotation of solubilized BR is discussed later. We further assume that deformation modes of helices, due to a simultaneous loss of some H-bonds, hardly occur within the EPR timescale (<200 ns). The high reestablishing capability of such hydrogen bonds should at least prevent distinct perturbations of the motional modes of R1. Hence, the position restraints applied to N, C, O were additionally used to properly stabilize the peptide planes and hydrogen bonds within the helical backbone. On the other hand, the side chains and all backbone atoms of the core of the EF loop (residues 160–164, Fig. 1 a) were kept unrestricted to allow free conformational sampling. RMSD analyses revealed that the nonrestriction of -carbons of helical amino acids do not significantly perturb the backbone structure. Further, the flexible -carbons are expected to distinctly improve the authenticity of side-chain motions.

View larger version (60K): [in this window] [in a new window]   FIGURE 3  (a) Front and side view of the model system containing a cysteine-bound MTS spin-label (black sticks) within a helix of 15 glycine amino acids (black ribbon). The contour surfaces cover the space that can be occupied by the nitroxyl nitrogen atom NX in the headgroup of the spin-label during a 6-ns, in-vacuo MD simulation at 600 K. The enclosing transparent area (light shaded, smoothed) limits the entire accessible space of the nitrogen, whereas the inner contour (dark-shaded, smoothed) qualitatively indicates the area with distinctly higher residence probabilities, which are located close to the peptide surface. (b) Each line or solid graph represents the accumulated covered volume of NX during different MD simulations. The covered volume was calculated as the sum of voxels (size: 1 Å3) sampled at least once during the prior simulation time. Results of MD simulations (600 K, in vacuo) of the spin-label bound to the most exposed positions within the E-F loop of BR (161–165) are given in the left panel. In increasing order of the final value: 163 (black solid), 164 (dark-shaded solid), 161 (shaded solid), 165 (light-shaded solid), and 162 (black line). The right panel shows the covered volumes by NX within the model system obtained under three different simulation conditions: 600 K, in vacuo, 20 ns (shaded line); 300 K, in vacuo, 40 ns (black line); and 300 K, explicit water, 50 ns (light-shaded solid graph). Both 300 K plots are sums of two independent simulations with -S-S- dihedral angles 90° (first half) and 270° (second half) each. The vertical lines at 20 ns and 25 ns mark the transitions between these simulations. The ordinates of both figures are identical. Specific conditions of the 600 K MD simulations are given in the text.

Apart from the final spectra calculation the mean-square fluctuation of the Euler angle ß, mainly reflecting the reorientation of the nitroxide -orbital, has been used to estimate the mobility of R1 (22). A more detailed comparison of simulated and experimental mobility parameters are given in Results and Discussion.

The potential energy topology
An effective potential energy function that determines the nitroxide reorientational motion was calculated from the simulated orientation probabilities. For this purpose the Euler angle space was divided into 500,000 unit cells () enabling a resolution of 3.6° in , ß, and . This resolution has been chosen to assure that the orientation distribution of the most rigid spin-label position led to the occupation of at least five cells in each angle dimension. The hits within every unit cell () were counted, normalized with respect to the total number of hits, and assigned to the center of the cell. The three-dimensional histogram, thus obtained, reflects the complex population distribution of the nitroxide orientations. Due to the limited sampling time of 6 ns, the topology is distinctly modulated by noise. This noise can cause large fluctuations of the potential-dependent motion during the single particle simulation and should be eliminated without influencing the global topology. Smoothing algorithms that utilize averaging techniques tend to an undesirable broadening of sharp peaks caused by that noise. To avoid this problem we applied the MEDIAN smoothing algorithm (23) with the smallest filter matrix (3x3x3) to fully eliminate the sharp noise peaks. The perturbation of the global topology due to the MEDIAN filtering has been proven to be less than that caused by an averaging square filter. Since the MEDIAN technique can only substitute original values with already existing adjacent values, it tends to form small steplike plateaus. However, this effect is less significant when smoothing is performed in high dimensions. If necessary, an additional square filtering could be applied to straighten out these remaining plateaus. The edges of the grid can be easily smoothed taking advantage of the periodic boundary conditions in , ß, and .

In the next step, the probability topology W() was converted into a free energy topology E_eff(),

(1)

where () denotes the , ß, triple of the centers of the unit cells; k_B is the Boltzmann constant, T the temperature (300 K), and C an arbitrary constant.

The single particle simulation
The calculation of a proper EPR spectrum requires reorientation trajectories with lengths of a few hundred nanoseconds (8). Such trajectories were obtained by the numerical solution of the potential-dependent Langevin equation for a single particle undergoing rotational diffusion (24):

(2)

The equation provides three contributions to the general torque M_i of a particle having a moment of inertia I_i. The first term on the right side represents the friction torque on a particle with the drag coefficient _i and an angular velocity _i(t) in the x, y, or z direction. The drag coefficient is related to the rotational diffusion coefficient D_i by the Einstein equation:

(3)

The expression T_i((t)) is the torque due to an effecting potential and is given by the gradient of the free energy topology for an orientation :

(4)

The value of E_eff at a distinct orientation is obtained from the discrete free energy topology E_eff() by trilinear interpolation in the Euler angle space. By means of the last term R_i of Eq. 2, a stochastic torque is implemented to mimic the coupling with the solvent. The value R_i satisfies the fluctuation-dissipation theorem

(5)

where _ij is the Kronecker and (t) the Dirac delta function. According to the proposed algorithm (20,25) for solving Eq. 2 and for I_i/_i << t, the reorientation _i(t_n+t) is given by

(6)

with the torques T_n,i = T_i( (t_n)). The value T_n,i is further assumed to be approximately constant during each time step. Thus, the derivative of T_n,i can be neglected. Under consideration of Eqs. 3 and 4 we get a conventional Brownian dynamics equation:

(7)

The last term _n,i of Eqs. 6 and 7 is a random variable that reflects the contribution of R_i. The value _n,i meets the rules of a random walk with

(8)

The probability distribution P_i of _i is given by the Fick's second law leading to a Gaussian distribution with zero mean and a mean-square value of 2D_it:

(9)

The arithmetic implementation of this Gaussian probability distribution was made using a normally distributed random number generator.

We generated a set of 160 single particle trajectories for each label position with lengths of 720 ns (18,000 data points per trajectory). The time step per iteration of the simulations was 2.5 ps. Starting orientations were obtained from randomly chosen angle triplets of the Euler angle trajectory of the MD simulation. The average W() histogram of these 720-ns trajectories showed an excellent agreement with that of the original MD trajectory for all analyzed R1 positions.

Spectra calculation
The calculation of the EPR spectrum is based on the approach described in detail by Steinhoff and Hubbell (9). Due to the anisotropy of the g- and A-tensors, the Larmor frequency depends on the orientations (t) of the nitroxide group with respect to the external magnetic field B₀, which is defined to be parallel to the z-direction of the laboratory frame. The Larmor frequency trajectory is determined from the energy eigenvalues of the spin Hamiltonian in the laboratory frame. An appropriate expression of the spin Hamiltonian is

(10)

with g-tensor g and hyperfine tensor A, both diagonal in the nitroxide frame. The value ß_e is the Bohr magneton, S and I are the electron and nuclear spin operators. This Hamiltonian is converted from the nitroxide frame into the molecular frame of the protein using a transformation matrix L, which is a function of the Euler angles , where index t denotes transposition:

(11)

Considering the protein reorientation with respect to the B-field, a second transformation matrix L' is necessary to finally transform the Hamiltonian into the laboratory frame. In the case of BR embedded in membrane sheets, the rotational correlation time _BR of the protein is distinctly larger than the upper limit of the sensitive timescale of X-band EPR of 200 ns. Hence, the molecular frame can be treated as static within the laboratory frame in the presented EPR spectra simulations.

To determine the eigenvalues of the Hamiltonian, two convenient approximations of Eq. 11 are utilized taking into account that the local nuclear magnetic field is small compared to the external field. The electron spin is, therefore, mainly quantitized along the z axis, whereas the x- and y-contributions of S can be neglected (intermediate field treatment). In the absence of motion, the Hamiltonian is then given by

(12)

with the expressions for g_zz and A_ij given in Steinhoff and Hubbell (9).

An approximation for fast reorientational motion assumes that the influence of the electron magnetic dipole on the nuclear spin becomes less significant in relation to the influence by the external magnetic field (26). Both the electron and nuclear spin are then quantized along the axis of the external magnetic field. The pseudosecular terms S_zI_x and S_zI_y of Eq. 12 vanish (high field treatment). Considering the autocorrelation function of the Euler angle ß, a weight function is applied to define a correct approximation between the two motional limits for a given spin-label mobility. This function facilitates a smooth transition between the high-field and intermediate-field treatments. A detailed description of this approach is presented elsewhere (9,26,27). From the energy eigenvalues, the energies of the allowed transitions are calculated for the spin-label orientations (t), which finally yield the Larmor frequency trajectory (t). For the presented spectra calculations, two different datasets for the principal values of g and A were applied, depending on the polarity of the predominant environment of R1 during the MD simulation. According to fits of several low-temperature EPR spectra of well-characterized positions in BR (e.g., 46, 100, 101, 129, 171, 193), we assigned principle values for water-exposed spin-label positions (g: 2.0083; 2.0065; 2.0026 and A: 0.56; 0.48; 3.65 [mT]) and for protein intrinsic R1 orientations (g: 2.0085; 2.0065; 2.0027 and A: 0.57; 0.41; 3.58 [mT]). A data set containing 320 trajectories of the complex magnetization M₊(t) was subsequently calculated from (t) by a numerical recursive, approximate solution of the Bloch equation:

(13)

This equation includes the assumption that the Larmor frequencies can be discretized over a small time step t, equal to the time step of the Euler angle trajectory (40 ps). The assumption was initially discussed by Itzkowitz (26) and indirectly justified via comparison of resulting EPR spectra with the results of Freed's MOMD model (9). To consider additional line broadening, the magnetization trajectories M₊(t) were multiplied by exp(–t/T₂). The relaxation time T₂ depends on the local environment of the spin-label. However, the variation is small, and T₂ values in the range of (20 ± 2) ns revealed agreement between simulated and measured spectra in almost every case. In addition, a convolution with a Gaussian function may be applied to account for inhomogeneous line broadening, e.g., due to unresolved hyperfine interactions with the methyl protons of the nitroxide ring.

To simulate the spectrum of a sample of randomly oriented spin-labeled protein molecules, the final magnetization trajectory was calculated as an average of 320 trajectories of a representative set of 256 different protein orientations (82,000 trajectories). The corresponding magnetization amplitudes were scaled by the sine of the angle between the protein z-axis and the magnetic field. The EPR absorption spectrum is then given by the real part of the Fourier transform of the averaged magnetization trajectory. To allow a convenient comparison of simulated and experimental spectra, all spectra are presented as the first derivative of the related EPR absorbance spectrum.

An EPR-related mobility parameter determined from MD simulations
A basic treatment is given here to obtain an EPR-related mobility parameter directly from the MD trajectories. Due to the definition of the nitroxide coordinate system (z-axis parallel to the p-orbital of NX) and the rotation sequence of the Euler transformation (around z, y', z''), the Euler angle ß specifies the angle between the orientation of the p-orbital of NX and the magnetic field vector, fixed to the z-axis of the laboratory frame. To allow comparisons between MD simulations showing different mean orientations of R1, the z-axis of the laboratory frame is further arbitrarily defined to be parallel to the average orientation of the p-orbital of NX (ß' = ß–ß). The population distributions in the Euler angle space are then centered around ß' = 0 for each spin-labeled mutant. This is reasonable under the assumption that the distributions of molecule orientations are isotropic in all measured BR samples. The hyperfine tensor of the used spin-label is nearly axially symmetric and the principle values can be approximated by A_xx = A_yy = A and A_zz = A_||, where A_|| is 7–8 times larger than A for R1. Since we focus our attention to the spectral contributions of the largest spectral splitting, the spatially averaged value of A_|| in the laboratory frame has to be considered. According to Steinhoff and Karim (28), this value can be approximated for rapid reorientational motion by

(14)

As discussed by Steinhoff and Karim (28), Eq. 14 includes the simplification that the main orientations of the p-orbital of NX are roughly isotropically distributed around the magnetic field vector or show rotationally symmetric characteristics. This is treated to be reasonable because most of the investigated R1 sites revealed roughly spherical or ellipsoidal topologies of the mainly populated orientation states in the Euler angle space. The variance of the cosine of ß' (Eq. 14) is used here as an EPR-related mobility parameter of the MD simulations.

Solubilization of BR
The investigated samples of membrane-embedded BR mutants were prepared by M. Pfeiffer according to Pfeiffer et al. (12). Solubilization of selected mutants in Triton X-100 was performed as follows.

After 5 min of ultra-sonification to initially reduce the size of the membrane sheets, extraction of BR from the membrane-embedded state was performed by adding at least 300 Triton X-100 detergent molecules per BR monomer corresponding to a concentration of 6 mM Triton X-100. This procedure causes the loss of all trimer-trimer and monomer-monomer contacts of BR (29) and inhibits the ability of lipids to form large membrane bilayers. Additional EPR spectra of spin-labeled BR samples with higher excess of Triton X-100 (6–20 mM) were obtained to investigate the structural stability of solubilized BR. Triton-containing samples were prevented from intense radiation to avoid photobleaching. After the solubilization procedures, a characteristic shift of the absorbance maximum of the protein-bound retinal from 570 nm to 553 nm (30) was observed. Precipitation of solubilized BR mutants was not observed after prolonged centrifugation (1 h with 16,000 g).

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSION ACKNOWLEDGEMENTS REFERENCES

 
A representative selection of recently published structures of BR is superimposed in Fig. 1 b. The major part of the protein structure is well defined by x-ray and cryo-electron microscopy studies. Although the published BR structures agree well in the transmembrane region of the protein, including positions 165–171, they show distinct differences particularly in the conformation of the cytoplasmic loop between helices E and F. Both presented cryo-EM models (Fig. 1 b; 2AT9, 1CWQ) represent the highest conformational diversity of the loop. Especially the loop structure presented in 1CWQ (31) reveals at least four positions with distinct discrepancies between the predicted degrees of reorientational freedom of R1 and the measured EPR spectra. The expected mobility of R1 bound to the two most exposed E-F loop positions (160, 161) of this structure disagree with the respective EPR spectra, which indicate relatively immobilized R1 side chains, whereas the partially buried positions 164 and 165 show EPR spectra of high R1 dynamics. Furthermore, the low collision frequency between the nitroxide attached to positions 157–161, 167 and the polar quencher chromium-oxalate do not correspond to the water-exposed residue orientations in this loop conformation (12). Most of the x-ray structures are characterized by uncertain E-F loop structures due to diffuse or absent electron densities in this region (13,32,33). However, the loop structure published by Essen et al. (1BRR) reveals no obvious discrepancies to the EPR spectra (13,34). This loop structure is furthermore almost identical to the E-F loop of the structures 1QHJ (35) and 1FBB (33) (Fig. 1 b) and the structure 1IW6 (36) (not shown). The present EPR spectra simulation algorithm is not expected to reveal a clear preference for only one of these very similar conformations. Hence, in this study a detailed comparison between experimental and simulated spectra on the basis of the conformation reported by Essen et al. (13) is performed to elucidate the capabilities of the presented EPR spectra simulation method. Since this structure consists of a complete trimer, the adjacent monomers were involved in the MD simulations as well. Approximately 10 C-terminal amino acids could not been adequately resolved in any model of BR, thus they were not considered.

MD simulations at different temperatures
The MD simulations were performed to sample solely the complex anisotropy of the reorientational motion of the nitroxide ring. It has been shown by LaConte et al. (37) that in vacuo MD simulations at 300 K are sufficient to calculate high order parameters, reflecting a high degree of immobilization of the spin-label. Although the applied simulations seem to yield realistic librational motion of the bound spin-label FDNASL, the question of whether the complete accessible conformational space is sampled within the 14-ns simulation time has not been addressed.

Here, we provide evidence that R1 side chains in exposed loop positions reveal much larger motional amplitudes in the time-range of EPR sensitivity. Special simulation parameters were chosen to focus on an entire conformational scanning, in which all significant restrictions are included. R1 bound to the loop position 165 of BR was selected to analyze the sampled Cartesian and Euler angle space at different temperatures, where the covered Euler angle space represents the reorientational freedom of the nitroxide. This initial investigation is based on the three-dimensional topologies of population densities calculated for each temperature (data not shown). Position 165 provides a solvent-exposed location and considerable opportunities for nitroxide contacts to the protein surface. Fifteen in-vacuo MD simulations (20 ns each) in the temperature range from 300 K to 900 K were performed, and the results were confirmed by simulations (300–750 K) for further positions with minor spatial restrictions (E161R1 and S162R1). We observed a drastic increase in the covered Euler angle space up to 550 K followed by only slight changes of the angle distributions between 550 K and 700 K. More frequent transitions over activation barriers of subsets of attractive van der Waals interactions between the nitroxide and the protein surface at 500–550 K, in-vacuo, are suggested as an explanation for the observed behavior. This enables a less hindered diffusion of the nitroxide headgroup along the protein surface. In fact, the revealed orientation distributions below 550 K are narrow and shift drastically their locations within the Euler angle space induced by temperature change. As already known, simulations at 300 K without explicit water molecules show a distinct tendency to enhance the internal density of the protein and exhibit an unrealistically strong surface tension. This effect is due to the presence of only protein-intrinsic, therefore anisotropic, nonbonded interactions and the absence of compensatory, mainly attractive protein-water interactions (e.g., H-bond and attractive van der Waals forces). Hence, particularly peripheral side chains, like those of the E-F loop amino acids, show an obviously restricted motion in 300 K in vacuo MD simulations. The observed restrictions seem to be too large even in relation to the hydrophobic effect resulting from real water environments. Increasing the temperature to 600 K weak short-range forces, like attractive van der Waals interactions, are assumed here to lose significance due to the related higher kinetic energy applied to each atom. In fact, on the basis of the used Lennard-Jones potentials, an increased mean separation of only 0.1 nm between the R1 headgroup and the protein surface could roughly halve the energy necessary to break a set of attractive van der Waals interactions. Above 550 K a global distribution is established with almost constant positions of the highest population densities for all analyzed R1 binding sites. In particular for in vacuo simulations, the degree of surface association of exposed side chains at 600 K is suggested to accomplish a more realistic behavior than at 300 K. The presented suggestions are substantiated by a detailed inspection of the side-chain dynamics at 300 and 600 K, given in the following sections. Above 700 K a tendency in the reorientational motion of R1 becomes dominant toward more uniform population densities in the entire accessible conformational space. At these temperatures, attractive van der Waals interactions are assumed to lose the ability to bias the motion of the spin-label.

Transition rates and EPR timescale
Although we neglect all kinetic information of the MD simulations for EPR spectra simulation purposes, an analysis of the transition rates is nevertheless important to estimate whether the populations of the found conformational states of R1 are equilibrated within the MD simulation time. Moreover, additional restrictions have to be implemented, if essential potential barriers are too small with respect to the kinetic energy at 600 K. The EPR spectral shape of a nitroxide at X-band frequencies is sensitive to reorientational motions with correlation times of 100 ps to 200 ns. Analyzing the reorientational dynamics of the R1 side chain during simulations at several temperatures, we identified two categories of transitions that are not within the EPR time window but nevertheless affect either the general diffusion rate or the maximum reorientation amplitude of R1.

The first category includes transitions over small energy barriers (in the order of 10 kJ/mol) established by individual attractive van der Waals, H-bond, and most of the Coulomb interactions involving the R1 side chain. An explicit-water MD simulation of R1 utilizing the OPLSAA force field and the TIP4P water implementation of GROMACS (38) was performed to study the lifetime of H-bonds between the nitroxide oxygen and surrounding water molecules. Using a polarization of the nitroxide group of OX: –0.3; NX: 0.08; C3/C4: 0.11 (partial charges) we could estimate an average lifetime of <1 ps. This value correlates with a transition barrier of 6 kJ/mol estimated by the Eyring equation and considering a symmetric bidirectional decay of the transition state. Due to the weak polarity of R1, the lifetimes of such individual interactions are usually shorter than 100 ps (barrier height of a transition between two states below 18 kJ/mol). Hence, we assume that a further increase in such high transition rates due to an MD simulation at 600 K do not alter the reorientation potential of R1 obtained by the MD simulation. However, the frequent appearance of these interactions leads to a significant contribution to the friction on R1 (e.g., diffusion along the protein surface). This effect is considered macroscopically using a realistic diffusion rate during the single particle simulation. Nevertheless, the cumulative strength of their collective occurrence, especially in the case of van der Waals interactions, causes the main restrictions even to the motional amplitude of the spin-label within the time range of EPR sensitivity. This will later be analyzed on a more macroscopic level and compared with the features of the hydrophobic effect. Strong Coulomb interactions could significantly bias the EPR spectrum. Hence, two different sets of partial charges of the nitroxide atoms (see Methods) were investigated empirically. Although surface-associated partial charges already appear too strong in in-vacuo MD simulations due to the neglecting of a surrounding dielectricum, the larger charge separation (OX: –0.3; NX: 0.08; C3/C4: 0.11) revealed slightly more consistent results for R1 in polar protein environments. However, further investigations might be necessary to optimize the settings. It is noteworthy that the tested charge separations have the same order of magnitude as the previously mentioned uncertainty in the quantum-mechanical determination of the partial charges for the nitroxide atoms.

The second category includes all strong restrictions concerning bond lengths, angle bending, repulsive van der Waals forces, etc. which cannot be overcome within 200 ns at 300 K (beyond the rigid limit). These barriers determine the principal restrictions of the R1 motion. Due to the very steep slope of such potential barriers, the gain in conformational freedom during MD simulations at 600 K with respect to 300 K is relatively small. Only slight corrections of the MD potentials were necessary to avoid unrealistic motion amplitudes at 600 K. The SHAKE algorithm was used to eliminate variations of all bond lengths. Improper dihedral potentials stabilized, e.g., the nitroxide headgroup and all amid planes, whereas weak position restraints, applied to the backbone atoms N, C, O of all helices, sustained the integrity of the helical H-bond network. Finally, a relatively high charge separation, especially in the nitroxide group, was used to facilitate sufficient Coulomb interactions involving R1.

To address the question of whether the helical backbone fluctuations significantly influence the spectral features, we determined the linewidth of the most buried positions on helix E (154,155) and helix F (170, 171). Experimental values of between eight and nine Gauss, which are similar to the linewidths of powder spectra, indicate a very rigid protein environment. Especially for position 171 the effective A_zz value of 3.47 mT measured at room temperature was found to be close to the A_zz value of the solid state (3.51 mT) obtained at 160 K. With the assumption that the backbone fluctuations through an extended part of a helix are of the same order of magnitude, the motion of the helical backbone in the vicinity of the E-F loop of BR can be considered negligible.

Intrinsic flexibility of the R1 side chain
An absolute prerequisite for a proper motion of the bound spin-label in 600 K MD simulations is a still-realistic intrinsic flexibility. To analyze the internal degrees of freedom of R1 under different conditions, we focus on the dihedral potentials—the weakest type of bonded interactions between covalently bound atoms. Fig. 2 presents histogram plots of the torsional states of all flexible bonds within the R1 side chain. The revealed torsional restrictions, based on the implemented dihedral potentials, are in addition partially amplified by, e.g., the angle-bending potentials and the 1–4 interactions. Using improper dihedrals, all atoms of the five-membered ring including the nitroxyl oxygen were restricted in the widely accepted coplanar conformation (39). MD simulations of our model system were used to compare the distributions of the essential dihedral angles under different conditions. The flexibility of R1 bound to a short glycine helix was determined from simulations at 300 K, both in vacuo (40 ns) and with consideration of 948 water molecules (50 ns), and the proposed in vacuo simulation at 600 K (6 ns). The results represent the dynamical behavior of an exposed protein-bound spin-label within a hydrophilic peptide environment.

In the case of the C-Cß (₁) bond, the familiar three-peaked plot was obtained showing the favored gauche/trans states of sp3-hybridized carbons. The peak positions are in perfect agreement with results of other ₁ dihedral analyses of native side chains, especially methionines and cysteines (40,41). A first indication for overestimated restrictions at 300 K in vacuo is that the 60° state hardly occurs during the entire 40-ns simulation time. Both -C-S- plots reveal almost ideally staggered conformations, which is also consistent with previously published results (42,43). The torsion angle distributions of the -S-S- bond show two characteristic states at 90° or 270°. In these cases the lone-pair electrons in the 3p atomic orbitals of both sulfur atoms are orthogonal to each other and their repulsion is minimized. An exhaustive survey about the electronic properties of the disulfide conformations is given by Boyd (44). As discussed later, transitions between these two low-energy orientations are unlikely to occur at room temperature within a few nanoseconds, whereas they were observed in all 600 K MD simulations (average lifetime is 1 ns). Hence, all presented histograms obtained at 300 K are 1:1 superpositions of two simulations with different initial values of the disulfide dihedral angle (90° or 270°). Apart from the high transition barrier, Jiao et al. (45) and Maxfield et al. (43) found that the energy penalty for distorting both ground states within 30° in each direction is small. The angle distributions determined here from 600 K MD simulations are in adequate agreement with these results. Rotation around the last bond -C-C1- shows the highest flexibility and most frequent transitions. With exception of the syn- and anti-periplanar conformations, almost all orientations are populated. If the entire R1 side chain is located in a more restrictive environment, the rotation around -C-C1- is dominant. Since this bond is nearly parallel to the nitroxide bond, it mainly affects the reorientational motion of R1 around the x-axis of the spin-label system.

Comparing the torsional distributions obtained under the three different conditions (Fig. 2), it is evident that the general positions of highest and lowest populated states are in good agreement. Even at 600 K the relatively weak dihedral potentials lead to acceptable anisotropic distributions. Nevertheless, as visible particularly for both -C-S- bonds, the sparsely populated energy-rich orientations are slightly more frequently occupied at 600 K. This indication of an increased internal flexibility of R1 does not necessarily result in a significantly higher orientational freedom of the nitroxide headgroup since particular orientations of the nitroxide can generally be achieved with different sets of dihedral angles. The dihedral angle distributions of R1 attached to exposed positions of BR display shapes similar to that of the presented 600 K MD simulation of the model system. However, the population ratios of the frequently occupied states differ due to the specific environmental restrictions. Buried or motionally highly restricted R1 side chains provide very narrow distributions of the dihedrals. A predominant steric hindrance can induce several torsion angles to occupy energy-rich values. Such limited gains in the potential energy of R1 are assumed not to be in contradiction to the behavior under experimental conditions.

Due to the implemented additional restrictions in the 600 K MD simulations mentioned previously, we assume that only a few conformational changes might occur during a 600 K simulation which are probably not observable by X-band EPR spectroscopy (lifetimes above 100 ns at 300 K). Referring to the time range of EPR sensitivity, the dihedral reorientation of the disulfide bond is the most critical remaining transition. Both energetic barriers between the two -S-S- ground states at 90° and 270° seem to be very high compared to other single bond dihedral potentials in proteins. Most of the published values are in the order of 25–40 kJ/mol and 37–50 kJ/mol for transitions via the trans and the cis conformations, respectively (44–46). For the preferred trans, transition lifetimes of several hundred nanoseconds can be estimated for both ground states at 300 K. On the other hand, the average lifetime determined by in vacuo MD simulations at 600 K amounts to 1 ns only. This is the result of seven trajectories of the most exposed label positions within the loop of BR. The -S-S- dihedral potential is approximated by the GROMOS force-field using a sinusoidal function with a potential barrier of 33.4 kJ/mol. Additional contributions are repulsion forces due to angle bending and the 1–4 Lennard-Jones interaction between Cß and C, which lead to an asymmetric potential with a less favored cis transition. Thus, the cis/trans energy barriers used by GROMOS seem to be close to the upper limits of the published values. The expected discrepancy in the dihedral flexibility of disulfide bonds at 300 K and 600 K cannot be fixed without distinct expense. A higher potential barrier may avoid transitions during the 600 K simulations. In this case we have to perform two independent simulations—one for each dihedral ground state. In addition, a third 600 K MD simulation with unaltered or lower barrier height would be necessary to obtain the fractions of R1 side-chain orientations at 90° and 270° within a macroscopic EPR sample. However, torsional transitions around the disulfide bond stimulate compensating collective reorientations of almost all other dihedral angles within the R1 tail. This is necessary to prevent drastic and sudden displacement of the R1 headgroup. Furthermore, a comparison of the covered space of the nitroxide nitrogen during simulations with pure 90° or 270° dihedral orientations shows large overlap. Hence, the total influence on the location of the nitroxide group is relatively weak; see also Fig. 4 (300 K in water simulations) and Fig. 5. According to these results we applied 600 K MD simulations without amplification of the -S-S- dihedral potential. The lifetimes of all other dihedral orientations, shown in Fig. 2, are distinctly shorter due to much smaller potential barriers.

View larger version (53K): [in this window] [in a new window]   FIGURE 4  The three parameters of the DTL system are illustrated on the structure of R1: The first parameter is the length of the vector CB-NX. The spreading angle is defined by the vectors CA-CB and CB-NX. The dihedral angle is the angle between the plane through the vectors N-C and CA-CB and the plane through the vectors CA-CB and CB-NX. , -histogram plots (TL plots) are efficient for revealing the reorientational behavior of a side chain. Here, TL plots (interval resolution: 1.8°) obtained from R1 of the model system are compared. The outer open areas in the TL plots were not sampled, whereas the gray-scale gradients indicate the histogram counts from black (at least one hit) to white (maximum amount of hits). The upper graph presents R1 orientations obtained from the 600 K, in-vacuo MD simulation (6 ns). Data from 300 K MD simulations with explicit water are visualized in the second line, whereas the results of 300 K, in-vacuo MD simulations are shown in the third line. For both 300 K conditions independent MD simulations for each dihedral ground state of the disulfide bond (90° and 270°) were performed (left and middle panel) followed by a superposition of both (right panel). Due to observed transitions between the ground states of the disulfide dihedral within the 600 K MD simulation (see Fig. 2 for the population ratio) the upper graph should be compared with the superpositioned (rightmost) panels of the 300 K MD simulations.

View larger version (15K): [in this window] [in a new window]   FIGURE 5  Simulated and experimental EPR spectra of the exposed R1 side chain bound to position 165 of BR. The simulated spectra in panels a–c were calculated based on the reorientational restrictions of R1 during long-time in-vacuo MD simulations at 300 K. Separate spectra simulations (20 ns each) were performed considering R1 conformations in one of both ground states of the disulfide dihedral, i.e., -S-S- of 90° (a) or 270° (b) (see text). Spectrum calculations based on the superposition of the orientation state populations of both 300 K MD simulations are presented in panel c. For each 300 K MD simulation (a–c), spectra are calculated considering three different rotational diffusion constants during the single particle simulations: D = 0.05 ns–1 (black dashed), 0.2 ns–1 (shaded line), and 0.5 ns–1 (black line). In panel d the measured spectrum (black line) and the simulated spectrum based on a 600 K MD simulation (shaded line) (D = 0.048 ns–1) are given for comparison.

Accessible sites for spin-labeling are usually in close vicinity to the protein surface. Apart from a small fraction of completely buried label orientations, the motion of the bound spin-label is therefore directly influenced by the hydrophobic effect. Less polar and mainly aprotic groups, like several side chains, perturb the hydrogen-bond network within the surrounding water layer. Thus, a driving force to reduce the surface/volume ratio leads to a compact aggregation of such groups. In particular, the distinctly apolar R1 side chain should generally be in close contact to the protein surface (47). The analysis of the sampled space for the nitroxyl nitrogen atom NX in the model system yields a qualitative representation of the mobility of R1 at 600 K (Fig. 3 a), which is in agreement with the expected motion behavior.

Due to a distinctly decreased diffusion rate at the protein surface in simulations with explicit water, an equilibrated orientation distribution had not been reached within 25 ns. In general, global modes of R1 motion are strongly perturbed by explicit water molecules. This is affirmed by the plot of the covered Cartesian volume of NX, which is a measure of the general sampling rate of R1 conformations, during both 25-ns MD simulations (Fig. 3 b, right panel). As already mentioned, transitions between both ground states of the disulfide dihedral were not observed in all 300 K MD simulations (maximum length: 25 ns). For positions without pronounced steric hindrance, R1 side chains with a disulfide dihedral angle of 90° or 270° should be almost equally represented in a macroscopic sample. Hence, simulations for both R1 species were done for every 300 K MD simulation. Since rare transitions of the disulfide dihedral are expected in long-time dynamics at 300 K, the covered volume of the consecutively summarized trajectories of both species is presented in Fig. 3 b (right panel). In contrast to the in vacuo MD simulations, a nearly constant level of the covered volume, indicating the shortest time to reach a conformational equilibration, was not achieved within the evaluated time window of 50 ns. Furthermore, a comparison of the sampled volumes during all simulations under different conditions clearly reveals that the final value of the covered volume at 300 K in water exceeds the highest level of the screened volume during the 300 K in vacuo simulations. A rough extrapolation of the maximum level for 300 K in water results in a value similar to those of the 600 K in vacuo simulation between 10 and 20 ns.

To perform a detailed inspection of the preferred conformational states within the covered space of R1, we introduce the DTL (distance, theta, lambda) coordinate system shown in Fig. 4. It is represented by

1. The length of the distance vector V between the atoms Cß and NX.
   
2. The angle between the C-Cß-direction and V.
   
3. A virtual dihedral angle , which is the angle between the plane defined by the vector N-C (backbone orientation) and C-Cß and the plane defined by C-Cß and V, respectively.
   

Using these linearly independent DTL parameters, the side-chain motion can be expressed in relation to the C-Cß bond and to the local backbone orientation. The length of V is generally in the range of 0.4–0.8 nm. Its value depends strongly on the torsional states of the -C-S- bonds. Histograms of the angle and the dihedral angle (TL plots) are sufficient to observe all possible orientations of R1 (Fig. 4). The analyzed TL plots of R1 of the model system reveal a high probability of NX positions located within a sinusoidal belt that is defined by the largest possible values of angle . This upper limit of is given by the steric hindrance due to the helical backbone. Therefore, the frequently occupied belt indicates R1 orientations with close surface contact of the R1 headgroup. It corresponds to the dark-shaded ring of mainly occupied NX positions in the Cartesian space (Fig. 3 a). The sinusoidal shape in the TL plots results from the fact that the plane of this circular distribution of mainly occupied NX positions is not completely perpendicular to the C-Cß bond. R1 orientations without contact to the helix are represented by angle values below the belt.

The two TL plots of the 300 K in-vacuo simulations (Fig. 4 below) illustrate the previously suggested strong motional restriction due to an exaggerated surface tension. Here, the simulation with a dihedral angle of 270° revealed the highest observed flexibility of R1 at 300 K in vacuo. Although the spin-label is mainly oriented in one direction, diffusion along the surface and weak sampling of most of the feasible orientations are nevertheless visible for both simulations. All performed 300 K MD simulations of R1 on more complex protein surfaces showed distributions similar to the 90° dihedral state behavior (Fig. 4, below-left panel) due to a much larger amount of buried conformational states with higher numbers of van der Waals interactions involving R1 than in the glycine helix model (data not shown).

The motion behavior during the MD simulations with explicit water showed a more uniform screening along the surface, indicative for significantly weaker interactions of R1 to local spots on the peptide surface. Further, the sampling rate of orientations apart from the surface was higher than that of the 300 K in vacuo simulations. In addition to the results of the covered volume analyses, another clear indication for a still not equilibrated conformational screening in the presented explicit water MD simulations is given by distinctly differing TL plots obtained after 16, 19, 22, and 25 ns simulation time. On the other hand, all in vacuo simulations show negligible changes in the TL plots during the last 25% of simulation time (data not shown).

As already discussed, transitions of the disulfide dihedral states have been observed in all 600 K MD simulations but not in 300 K simulations within 25 ns. Hence, superpositions of the obtained probability topologies of both dihedral states were considered for all 300 K simulations (Fig. 4, rightmost panels). These superpositions imply frequent transitions between 90° and 270° within the timescale of EPR sensitivity (200 ns). Since this behavior is not realistic, the real conformational restriction of R1 is in between the discrete and superpositioned topologies. Nevertheless, the TL plot obtained from the 600 K MD simulation should be compared with the superpositions (Fig. 4, rightmost panels). Even at 600 K the lack of compensational interactions, due to the absence of solvent molecules, generates a remaining force toward the surface of the peptide. Therefore, as in the presence of a hydrophobic effect, the major orientations of R1 are those with distinct surface association. On the other hand, the high kinetic energy counteracts unrealistic conditions of surface tension and solute density. The ratio of NX positions aggregated and separated from the peptide surface is smaller at 600 K compared to the 300 K MD simulations. However, we assume a similar ratio in the case of the 300 K MD simulations with explicit water when a complete conformational equilibration is achieved (see Fig. 3 b).

The main peak in the TL plot of the 300 K in-vacuo MD simulation (-S-S- = 90°) is identical to the position with the highest occupation probability of the 600 K simulation, but this position was weakly sampled in the explicit water MD simulations (Fig. 4). The area of the peak is characterized by the maximum possible values of and a large contact surface to the helix atoms (represented by the surrounding open area). An R1 side chain with such a conformation is wrapped round the glycine helix along the helical groove, yielding the highest number of van der Waals interactions between R1 and the atoms of the helix. The corresponding Cß-NX distances (first DTL parameter) show the highest variation (0.44–0.76 nm) in relation to regions with other -values. Hence, the mentioned area in the TL plots represents a large variety of surface-attached R1 conformations. Interestingly, a 600 K in-vacuo MD simulation with neglected Coulomb interactions (all partial charges of R1 set to zero) reveals a shift of the most probable R1 orientation away from the mentioned area to smaller -values. Thus, it seems that a combination of all nonbonded interactions possibly leads to overestimated attractive forces for R1 orientations with intense surface contact in in-vacuo MD simulations even at 600 K.

Due to the relatively high polarizability of sulfur atoms, van der Waals interactions under participation of the disulfide sulfurs are relatively strong. This fact is implemented in the MD force field via the Lennard-Jones parameters for sulfur. In addition, a relevant hydrogen-bond interaction between the S of R1 and the hydrogen atom bound to C has been suggested (48). A common method to estimate the strength of hydrogen bonds is a comparison of the measured atom-atom distances in low-temperature crystal structures with the corresponding sum of their van der Waals radii. If sulfur atoms are involved, this method might fail (49). Atoms, like sulfur, with occupied high level electronic orbitals (due to hybridization) and/or lone pairs are not well defined by an isotropic electron density distribution provided by the spherical van der Waals model. Thus, the Lennard-Jones potential should be considered as rather n