| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
* Laboratory of Molecular Biophysics, Department of Biochemistry, University of Oxford, Oxford OX1 3QU, United Kingdom; and Biomolecular Modelling Laboratory, Imperial Cancer Research Fund, London WC2A 3PX, United Kingdom
Correspondence: Address reprint requests to Mark S. P. Sansom, Laboratory of Molecular Biophysics, The Rex Richards Building, Dept. of Biochemistry, University of Oxford, South Parks Rd., Oxford OX1 3QU, UK. Tel.: +44-1865-275371; Fax: +44-1865-275182; E-mail: mark{at}biop.ox.ac.uk.
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSION AND CONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSION AND CONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
). Due to the high expression levels of pore-forming proteins and transporters that enable the diffusion of nutrients into the periplasm, the outer membrane is considerably more permeable than the cytoplasmic or inner membrane, so much that a permanent electrochemical potential cannot be maintained with respect to the cell exterior. Thus those processes mediated by the outer membrane that require an energetic supply, such as active transport or secretion, must be necessarily coupled to the inner membrane, where electrochemical gradients and/or ATP hydrolysis can be used for generating this supply. Uptake of iron siderophores and vitamin B12 into the periplasm has been shown to be an example of an energy-dependent translocation process through the outer membrane. Numerous mutational and biochemical studies (as reviewed in Kadner (1990), Postle (1990), Moeck and Coulton (1998), Braun and Braun (2002), and Faraldo-Gómez and Sansom (2003)) suggest that across a wide range of these systems involve the energetic coupling of high-affinity receptors residing in the outer membrane, such as FhuA and BtuB, with a protein complex (formed by the proteins TonB, ExbB, and ExbD) anchored in the plasma membrane. Although little is known about the actual mechanism of function of these systems, it is widely accepted that the TonB protein utilizes the proton motive force across the cytoplasmic membrane to induce the translocation of ligands from the exterior of the cell into the periplasm, by causing a conformational change in the outer membrane receptor, via a physical linkage with its periplasmic face; the extent of this structural change, however, remains very much a matter of debate.
From a structural viewpoint, the best-known TonB-dependent transporters are FhuA, FecA, and FepA from Escherichia coli; respectively, these proteins mediate the uptake of ferrichrome, ferric enterobactin, and ferric citrate. Although the atomic structure of FepA has been solved only in the absence of its ligand (Buchanan et al., 1999), the structures of FhuA and FecA have been determined in both the siderophore-free and bound states (Ferguson et al., 1998, 2002; Locher et al., 1998;), revealing allosteric changes thought to promote the interaction of the siderophore-bound receptor with the transducer TonB. Interestingly, the FhuA structure has also been determined in association with two siderophore-antibiotic conjugates (Ferguson et al., 2000, 2001), illustrating the utilization of the TonB-dependent uptake pathway for delivery of antibiotic agents into the bacterium.
The structures of FecA, FepA, and FhuA are remarkably similar, which suggests that a homologous topology might characterize other members of the TonB-dependent receptor family. They consist of large membrane-spanning monomeric ß-barrel domains, composed of 22 ß-strands of variable length, many of which extend significantly beyond the membrane hydrophobic core into the extracellular space. Adjacent strands are connected generally by short periplasmic turns and long extracellular loops. In contrast with porins, in FhuA these loops do not fold back into the barrel, but rather project away from the membrane surface. Instead, the barrel pore is occluded by a large, globular amino-terminal domain, termed the plug (Locher et al., 1998).
The plug domain effectively divides the protein into two regions: at the extracellular side, loops in both the plug and barrel domain form the ligand binding-site, in which the siderophores are oriented so as to facilitate the formation of hydrogen-bonds between its iron-chelating groups and protein residues; at the intracellular face, the N-terminus of the plug projects toward the periplasmic space. Interestingly, the extent to which this region extends into the periplasm is increased in the case of the ligand-loaded structures, due to the unwinding of the so-called switch helix, in response to what has been proposed to be the propagation of conformational changes along the plug domain (Ferguson et al., 1998, 2002; Locher et al., 1998; Merianos et al., 2000; Coggshall et al., 2001). It has been suggested that this conformational change at the N-terminal domain of the outer membrane receptor might be required for TonB to interact preferentially with ligand-bound receptors, given the larger proportion of TonB-dependent transporters with respect to TonB complexes (Braun, 1998).
Molecular dynamics (MD) simulations of a membrane protein embedded in a fully solvated phospholipid bilayer have been used to study a variety of membrane proteins ranging from single transmembrane (TM) 
-helices to complex ß-barrels and -helix bundles (Forrest and Sansom, 2000). Recent simulations include those of bacterial outer membrane proteins such as OmpF (Tieleman and Berendsen, 1998; Im and Roux, 2002) and OmpA (Bond et al., 2002), of bacterial ion channels including KcsA (Guidoni et al., 1999, 2000; Bernèche and Roux, 2000, 2001; Shrivastava and Sansom, 2000, 2002; Shrivastava et al., 2002) and MscL (Elmore and Dougherty, 2001; Gullingsgrud et al., 2001), and of members of the aquaporin family of transmembrane pores (de Groot and Grubmüller, 2001; Jensen et al., 2001; Tajkhorshid et al., 2002).
We present here the results of MD simulations of duration 10 ns of the FhuA protein embedded in a hydrated membrane-mimetic environment (dimyristoyl-phosphatidylcholine lipid), in both the ligand-free (FHUA, Protein Data Bank (pdb) code 1by3) and bound states (FHUA+F, pdb code 1by5) (Locher et al., 1998). In particular, we will start by describing the similarities and differences between our model systems and the crystal structures, in terms of structural drift and secondary structure. Our simulations will be further characterized by analyzing the interactions observed between the plug and barrel domains, as well as between the siderophore (ferrichrome) and the protein. Next we will report on the most relevant differences observed between the two simulations, namely in the motion of loop L8 at the extracellular side, and in the permeation of water across the barrel interior; previously reported allosteric changes upon ligand-binding will also be discussed. To conclude we will analyze the atomic motions in the simulations in terms of calculated Debye-Waller factors and will discuss the extent of sampling of these motions during the simulations.
|
METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSION AND CONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
) molecular dynamics simulation package (www.gromacs.org—also see Lindahl et al. (2001)). The systems to be simulated were set up using a rectangular box (see Fig. 1 A) employing periodic boundary conditions. Numerical integration of the equations of motion used a time step of 2 fs, with atomic coordinates saved every 10 ps for analysis. Both simulations were conducted at a constant temperature of 310 K and a constant pressure of 1 bar (Berendsen et al., 1984). Solvent (i.e., 100 mM salt solution plus counterions), lipid, and protein were coupled separately to the temperature bath, with a coupling constant of 
T = 0.1 ps. Anisotropic pressure coupling (i.e., independent in the x, y, and z directions) was used, with a coupling constant of P = 1.0 ps and a compressibility of 4.5 10-5 bar-1. Energy minimizations were performed using a steepest descent algorithm. When restraints were required on, e.g., protein atoms, a harmonic potential with a force constant of 10 kJ mol-1 Å-2 was applied.
|
|
|
|
) was used, as implemented in GROMACS; aliphatic groups in both the protein and lipid were treated using the united atom representation. Bond distances were constrained using the LINCS algorithm (Hess et al., 1997). The van der Waals interactions were modeled using a 6–12 Lennard-Jones potential, cut off at 10 Å. Electrostatic interactions were computed using a particle mesh Ewald algorithm (Darden et al., 1993), with a 9 Å cutoff for the direct space calculation; the reciprocal space calculation was performed using a fast Fourier transform algorithm. The calculation of nonbonded interactions used first- and second-neighbor exclusions and scaled-down third-neighbor coefficients. Both simulations were set up on an O2 SGI workstation (Silicon Graphics, Mountain View, CA) and performed on an SGI Origin 2000 (Silicon Graphics), using eight parallel 195 MHz R10000 processors; the CPU time was 
150 h per ns of simulation.
Construction and equilibration of the model systems
A lipid bilayer was used to mimic the membrane environment in which FhuA exists; this lipid bilayer, initially made up of 288 dimyristoyl-phosphatidylcholine (DMPC) lipids, had been equilibrated during 1.5 ns of MD, resulting in an area per lipid of 62.4 Å2 (cf. experimental values of 60.0 Å2 (Petrache et al., 2000) and 59.7 Å2 (Petrache et al., 1998) obtained at 303 K by 2H NMR and x-ray diffraction, respectively). A 100 mM NaCl solution (SPC water model) was also used, alongside a number of counterions to neutralize the total charge of the system (all ions were placed initially at random). To set up the protein-membrane systems, a suitable cavity in the interior of the bilayer was generated using the protein solvent-accessible surface as a template; this methodology has been described in more detail in Faraldo-Gómez et al. (2002). A topology for the siderophore ferrichrome was generated ad hoc for the simulations of FHUA+F. Parameters were adapted from molecular mechanical calculations by Zinelabidine et al. (1993), and from GROMOS87 (Hermans et al., 1984) and QUANTA/CHARMM (www.accelerys.com) force fields. A 1 ns simulation of the free molecule in SPC water was used to validate this topology. Calculations of the pKA of ionizable residues were performed using the University of Houston Brownian dynamics solver of the Poisson-Boltzmann equation (Davis et al., 1991) as well as in-house code; this procedure has been described in more detail in Bashford and Karplus (1990) and Adcock et al. (1998).
The protein-bilayer systems (71,000 atoms) were equilibrated for 1.5 ns of MD, so as to improve the packing of lipids around the protein and the distribution of the ions in solution; the resulting C root mean-square deviation (RMSD) of the transmembrane domains, with respect to the crystal structures (1by3.pdb and 1by5.pdb, (Locher et al., 1998)), was 0.9 Å in both cases. The resulting protein structures were then replaced by the crystallographic conformations and the system energy minimized, after which 50 ps MD runs with positional restraints on the protein atoms (except hydrogen atoms) were conducted. The production runs thereafter were of 10 ns duration each.
Analysis
The secondary structure of the protein snapshots during the simulation was analyzed with the DSSP algorithm (Kabsch and Sander, 1983); the results are shown in terms of frequency of occurrence of either ß-strand or -helical conformation. Standard and adapted GROMACS routines as well as in-house code were utilized in the analysis of RMSD, mean-square fluctuations (MSF), hydrogen bonding, interatomic distances, and water residence times. When analyzing structural drift with RMSD calculations, for instance in loop n in the barrel domain, "self-fit" refers to a calculation in which the RMSD of loop n is calculated after least-square fitting the protein snapshots using loop n as the reference, thus removing the rigid-body motion of loop n from the analysis. In contrast "TM-fit" refers to a calculation where the snapshots are fitted using the transmembrane region of the ß-barrel as the reference domain, and is therefore useful for detecting swinging motions.
Conformational sampling was assessed using a block analysis of the MSF of the C atoms. The MSF for a certain subset of C atoms and a time window w is defined as:
 |
Xi 
is the average of the latter over all the snapshots in the window. Thus, the window-averaged MSF for the subset is simply:
 |
and a simulation time T is W = T-1. Finally, simulation B-factors are deduced from the expression B = (8
2/3) MSF, and therefore carry the same dependence on the sampling window length as the MSF. |
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSION AND CONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
|
|
|
-helical segments are seen in the loop L4, and to a lesser extent in loop L3, which are the longest loops in the barrel domain. Overall, the simulation secondary structure profiles closely resemble those of the crystal conformations. Significant changes were observed in localized regions (e.g., ends of ß8, ß9, ß13, and ß14, as well as in the small helix in L3), although these appear to result in a higher degree of convergence in the secondary structure profiles of FHUA and FHUA+F, compared to the initial conformations.
Analogous conclusions can be drawn for the plug domain (Fig. 2, C and D), i.e., that the most significant secondary structure features remain stable throughout both simulations. This domain contains a well-defined four-stranded ß-sheet core (ßA–ßD), with a small fifth strand (adjacent to ßB) present permanently only in the ligand-unloaded state. Stable -helical stretches are also seen in both simulations in the regions defined here as S3, S4, and S5, as well as at the C-terminal end of loop LC. Partially stable ß-hairpin-like motifs delimit the S2 region and loop LA. As can be seen, the most apparent difference between the ligand-free and loaded crystal structures of FhuA, that is, the fold of the so-called switch-helix in region S1, is maintained in the simulations; moreover, in FHUA this helix is somewhat lengthened at its C-terminus, relative to the crystal structure.
Structural drift: the ß-barrel domain
It has become commonplace in protein MD simulations to analyze the RMSD of the backbone or C trace from the starting crystallographic structure as a simple means by which the stability of the simulations may be assessed. The results of such analysis for the current simulations (Fig. 3 A) reveal some general trends also seen in simulations of other membrane proteins (see below). After 10 ns of simulation, the structural drift of the transmembrane region of the backbone (or rather, the domain that spans the hydrophobic core of the membrane) has clearly stabilized with respect to the crystal conformation, with RMSD values of around 0.7 Å in FHUA+F and 1.0 Å in FHUA. Inclusion of the nontransmembrane segments of the ß-strands yields slightly higher and less stable RMSD profiles, in the range of 1.1–1.3 Å. Lastly, for the entire ß-barrel, that is, including the extracellular loops and periplasmic turns, the final deviation is 2 Å, but the continuous increase indicates that it does not represent a stabilized configuration.
To identify the source of the continuous drift in the barrel domain observed in both FHUA and FHUA+F simulations, in Fig. 3, B and C, we show the time-averaged deviations of the loops (L1–L11) and the turns (T1–T10). Clearly, the most significant deviations in the barrel domain occur on the extracellular side in loops L4–L8; in particular, loop L7 and especially loop L8 show large deviations that do not stabilize after 10 ns, as can be seen in the time series (data not shown). At the periplasmic face of the protein, the majority of the turns show a higher degree of mobility (larger standard deviations), but all the respective time-dependent profiles (not shown) are fairly stable at around 2 Å, with the exception of turns T7 and T10 in the FHUA+F simulation, which yield continuously drifting RMSD profiles of averages 3.0 Å and 3.5 Å, respectively.
It is noticeable that the largest deviations, seen in L7, L8, and T7, appear to correspond to conformational changes mainly involving a drift of the loop or turn with respect to the transmembrane domain, rather than a change in the conformation of the loop/turn per se. To illustrate this, in Fig. 3, B and C, we also show self-fit RMSD values (see Methods). These segments show structural changes within the loops/turns of 1 Å or less, with reasonably stable time-dependent profiles (not shown). Thus to a first approximation, the loops and turns may be viewed as internally stable structural elements that are able to move with respect to the transmembrane barrel on a nanosecond timescale, although we have not fully sampled these motions. The possible functional relevance of these motions will be discussed below.
It is useful to compare these structural drifts with those seen in some other simulations of membrane proteins embedded in lipid bilayers. For instance, simulations of -helical bundles, such as those of the potassium channel KcsA that used the low-resolution crystallographic structure, have shown overall backbone deviations in the range of 1.9–2.5 Å after 1–1.5 ns of simulation (Bernèche and Roux, 2000; Shrivastava and Sansom, 2000); a low value (1.6 Å) has been obtained in a subsequent study of the 2 Å-resolution structure after 2 ns of simulation (Domene and Sansom, 2003); and recent studies of the aquaglyceroporin GlpF have reported deviations between 1.2 Å and 1.5 Å after simulations of duration 1 and 2 ns, respectively (de Groot and Grubmüller, 2001; Jensen et al., 2001). ß-barrel pore-forming proteins appear to provide a lower bound in terms of protein drift: in an early study of the OmpF trimeric protein, RMSD values of roughly 2 Å were obtained for each monomer after a 1-ns simulation (Tieleman and Berendsen, 1998), but more recent work on this same porin has shown a deviation of 1.4 Å after 5 ns (Im and Roux, 2002), similar to that found in a recent simulation of the C-terminal 171-residue ß-barrel domain of OmpA (Bond et al., 2002). In some of the examples given above, the RMSD analysis was also performed for just the transmembrane secondary structure elements, resulting in deviations that are roughly half that of the whole molecule (e.g., 0.7 Å for OmpA, 0.8 Å for OmpF, and 1.0 Å for KcsA). This is consistent with the results from FhuA and suggests that large degrees of structural drift may be, in general, associated with the extramembranous loops that are exposed to an aqueous environment.
Structural drift: the plug domain
FhuA differs from, e.g., OmpF and OmpA in having an internal plug domain. RMSD analysis (Fig. 4) for this domain reveals that the structural drift is rather low and stabilized in its central ß-sheet core (which separates the extracellular and periplasmic faces), which deviates from the crystal conformation by <0.6 Å after 10 ns; when the solvent-exposed regions are included in the calculations, the resulting deviations increase in both FHUA and FHUA+F, but to a lesser extent in the latter simulation (1.6 and 1.2 Å, respectively).
|
1.6–2.0 Å in the case of FHUA, whereas in FHUA+F the motions in the structure of the unfolded switch-helix appear to be predominant.
Ferrichrome and the binding pocket
Extramembranous loops are crucial in regulating the activity of TonB-dependent receptors and other outer membrane proteins such as porins, by providing suitable sites for substrate recognition, selectivity, or binding (Klebba and Newton, 1998; Koebnik et al., 2000). Consistently with this functional role, these loops appear to be capable of adopting a rather stable conformation, by virtue of their association with neighboring loops or with the barrel domains. In FhuA, the binding site for ferrichrome is formed by loops LA, LB, and LC from the plug domain and loops L3, L4, L5, and L11 from the ß-barrel. During the simulations, we observe multiple interloop hydrogen-bonding interactions that may contribute to the stability of this arrangement. L5 appears to secure the conformation of loops L3 and L4 by folding around them and forming several hydrogen bonds that involve both backbone and side-chain atoms (for instance, in FHUA+F, an average of 7.7 and 6.5 H-bonds were detected between these segments, respectively). Analogously, loops L3 and L11 seem to be mutually stabilized by two permanent backbone H-bonds. In the plug domain, LA, LB, and LC are also engaged in similar H-bonding interactions (e.g., 3.1 H-bonds on average were found between LA and LB as well as between LB and LC during the FHUA+F simulation).
The specific characteristics of the siderophore binding site in TonB-coupled receptors are believed to vary according to the chemical nature of the corresponding ligand. In FhuA, the pocket is lined by numerous hydrophobic residues, which are predominantly aromatic. In Fig. 5 A, we show a simulation snapshot of ferrichrome with the side chains of its nearest neighbors. The siderophore is oriented such that the iron-chelating groups face the protein interior, whereas the cyclic peptide, or tail, remains solvent-exposed. The tail of the siderophore is significantly more flexible than the head, as analysis of the mean atomic fluctuations during the simulation reveals (averages of 0.3 ± 0.1 Å and 1.0 ± 0.5 Å for the head and tail backbone atoms, respectively). The location of the iron-chelating groups within the binding site is locked by the formation of long-lived hydrogen bonds between three of the six oxygen atoms that coordinate iron and residues Arg81 (plug) and Tyr244 (ß-barrel) (Fig. 5 B).
Signaling across the plug domain
The crystal structures of the ligand-free and bound states of FhuA (Ferguson et al., 1998; Locher et al., 1998) provided clues as to the currently upheld mechanism by which siderophore binding might be signaled across the outer membrane, triggering the active translocation of the iron complexes into the periplasm. Upon binding of ferrichrome, the formation of ligand-protein contacts induces a shift of 1–2 Å toward the siderophore in the two short helices (in the S3 and S4 segments) that follow the first strand of the plug's ß-sheet core (ßA); this subtle conformational change was proposed to result in the destabilization and unwinding of the N-terminal, periplasmic switch-helix (in the S1 segment), yielding an extended segment that swings toward the opposite barrel wall, displacing its terminal C atom by 17 Å. This extended conformation might in turn facilitate the interaction of the so-called TonB-box (residues 7–12, not resolved in the crystallographic structures) with the TonB-ExbBD energy transducer complex. In the current study, we aim to provide insights into the dynamic properties of the plug, and in particular to investigate whether the differences observed between the crystal conformations of FhuA persist when dynamic fluctuations occur.
As suggested by the RMSD analysis above, visualization of superimposed simulation snapshots of the plug domain, with and without the ferrichrome bound, reveals this domain is overall rather inflexible (Fig. 6), with the exception of loops LA, LB, and LC in FHUA and the N-terminal domain S1 in FHUA+F. The conformations of helices S3 and S4 were found to be well preserved throughout the simulations. In FHUA, the RMS deviations of these helices with respect to the central ß-sheet were as low as 0.8 ± 0.2 Å and 1.1 ± 0.1 Å, respectively, whereas their intrinsic structure (as revealed by self-fits) deviated by only 0.23 ± 0.05 Å and 0.66 ± 0.07 Å. Similarly, in FHUA+F, the deviations of these domains were 1.0 ± 0.2 Å and 1.0 ± 0.2 Å relative to the central ß-strands, and 0.29 ± 0.08 Å and 0.60 ± 0.08 Å with respect to themselves. In Fig. 6 B, we have superimposed the conformations of S3 and S4 from both simulations, to illustrate the two different positions of these regions relative to the ß-barrel. The RMSD values of S3 and S4 from one simulation to the other are 1.9 Å ± 0.5 Å and 1.1 ± 0.3 Å, respectively, whereas their intrinsic (self-fit) structures differ only in 0.4 ± 0.1 Å and 0.5 ± 0.2 Å.
In view of these analyses, it appears reasonable to conclude that despite substantial thermal fluctuations at 310 K, these segments maintain different orientations relative to the ß-barrel. S3 in particular exists in a distinct orientation dependent on the binding state of the protein. Nevertheless, the question remains as to how these differences mediate the unwinding of the switch-helix.
Finally, in Fig. 6 C, we show snapshots of the N-terminal switch-helix from both simulations to illustrate the preservation of the helical fold in FHUA and the absence of a coil-to-helix transition in FHUA+F. It is clear that the unfolded state of this helix in FHUA+F is associated with a displacement of T7 toward the center of the protein. In fact, the conformation of T7 in the ligand-loaded structure is sterically incompatible with a folded conformation of the switch-helix, as can be seen by superimposing the molecular surfaces of FhuA in both binding states (not shown). The possible relevance of this conformational difference will be discussed briefly in the next section.
A functional role of the ß-barrel extracellular loops
In contrast to the plug domain, the reported analysis of the crystallographic structures of FhuA indicated no significant conformational changes in the ß-barrel domain upon ligand binding. Thus, this domain was suggested to provide a suitable scaffold for the protein to carry out its activity, but did not seem to participate in the signaling or energy transduction mechanisms.
However, analysis of the FHUA and FHUA+F simulations indicates the barrel domain might have a more active role. As discussed above, extracellular loop L8 is the segment with the largest structural drift during our simulations. After 10 ns, the RMSD of L8 with respect to the transmembrane region of the barrel domain had reached 6 Å in FHUA and 8 Å in FHUA+F. Its internal conformation, however, remained roughly stable, indicating a swinging motion, with (self-fit) time-averaged RMSD values ranging from 0.96 Å to 1.02 Å within the same simulation. This is in contrast with the crystal structures in that no significant conformational change in loop L8 upon ferrichrome binding was reported, in either the structures solved by Locher et al. (1998) or by Ferguson et al. (1998), nor can a substantial difference be seen when the structures from the two groups are compared. However, this loop is one of the most disordered segments of the protein. For example, in the ligand-bound structure solved by Locher et al. (1998), the backbone B-factors assigned to residues 553–558 of L8 range from 87 to 103. Similar values were assigned to this segment in the structure by Ferguson et al. (1998), as well as in the respective ligand-free structures. Furthermore, analysis of the protein packing in the crystals obtained by both groups reveals that this region is not involved in crystallographic contacts, which presumably reflects its intrinsic flexibility.
Although loop L8 drifts in both FHUA and FHUA+F simulations, the direction of the motion differs. Specifically, L8 drifts toward the siderophore in FHUA+F, but away from the siderophore binding-site in FHUA (see Fig. 7). This may be quantified if we plot the inverse of the minimum distance between residues in loop L8 and the ferrichrome, time-averaged over 2 ns intervals (Fig. 7 A). From this, it is apparent that the displacement toward the ligand of L8 is especially marked in the case of residues Ser556, Phe557, and Phe558. The two Phe rings, which in the initial conformation point away from the binding pocket toward the solvent, at the end of the simulation are in close contact with the siderophore tail, thus enhancing the hydrophobic character of the binding pocket. Similarly, the backbone of Ser556 becomes close enough to the ligand so as to form a hydrogen bond with residue Fer4 of the ligand (see Fig. 5 B). Given the limited duration of these simulations, and hence the fact that these motions are only partially sampled (see below), it is difficult to assess conclusively whether the above-described motions indicate that loop L8 has a role in closing the binding pocket of FhuA (in a similar way to that suggested for FecA on the basis of its crystal structure (Ferguson et al., 2002)). Nevertheless, the results are highly suggestive and enable useful speculation. First, we suggest that L8 may have a functional role in gating the FhuA protein upon ferrichrome binding. Second, we note that turn T7 (see Fig. 6 C), whose conformation appears to be coupled to that of the switch-helix, is connected via ß15 to loop L8. This suggests that the motions of L8 might be related to the destabilization of the switch-helix, and thus, that the allosteric effect of ligand-binding could act not only through the plug domain but also by using one or more of the transmembrane ß-strands.
Association of the ß-barrel and plug domains
Two possible mechanisms have been suggested for the mechanism of translocation of the siderophore through FhuA and across the outer membrane (reviewed in, e.g., Sansom, 1999). In one mechanism, it is suggested that the interaction between FhuA and the Ton complex might destabilize the high-affinity binding site for ferrichrome and cause a rearrangement of the plug such that a water-filled channel is formed. The siderophore, now released, would passively diffuse along this pore into the periplasm, perhaps aided by several polar residues located on the barrel wall, which are conserved throughout the family of TonB-dependent receptors. In a second possible mechanism, the interaction with TonB is suggested to provide a mechanical force large enough to pull the plug domain out from the interior of the barrel, the siderophore then being released from the binding site directly into the periplasm.
Our simulations certainly support the view that the plug domain must undergo a substantial conformational change to allow passage of the siderophore. Mapping the time-averaged number of hydrogen bonds between plug and ß-barrel onto the surface of the plug domain (Fig. 8 A) reveals a significant number of stable contacts between these two domains. This is true for both simulations, the total number of H-bonding contacts being 67 ± 4 in FHUA+F and 65 ± 4 in FHUA. This could suggest that a substantial activation energy would be required to entirely remove the plug domain from the barrel interior, and some investigators (Ferguson et al., 1998) consider that it would be large enough to prevent removal. However, we note that comparable numbers of H-bonds can be found in some dissociable complexes, e.g., 58 at the interface between glutaminyl-tRNA synthetase and its tRNA (Nadassy et al., 1999). Therefore, dissociation of the two FhuA domains, or at least induction of a conformational change sufficient to allow the permeation of the ligand, could be facilitated by the presence of numerous water molecules solvating the interfacial region (Fig. 8 B) in both simulations: the total numbers of barrel-water-plug H-bond contacts were 56 ± 6 in FHUA+F and 55 ± 5 in FHUA. The presence of these water molecules might be important in helping to reduce the activation energy barrier associated with barrel-plug dissociation. Moreover, in the relative motion of heavily hydrogen-bonded moieties, the apparent activation energy can be very much less than required to break all the H-bonds, due to lubricating water (as mentioned) and also H-bond switching.
|
7 Å, would not be possible unless there was a significant rearrangement of the side chains lining the interface.
In Fig. 8 C, we show the two orientations where the largest distances between the barrel and the plug can be found. In the 60° view, in both the ligand-free and bound states, a large cavity with a maximum diameter of 8 Å can be seen underneath loop LA, as well as in the space between the N-terminal segments S1 and S2 and the ß-sheet core. However, LA itself and the segment between S3 and S4 would prevent access to these regions from the binding site. In the 240° view, analysis of the ligand-free plug reveals a pore-like cavity, lined by domains S2 and S5 on one side and the N-terminal half of LC on the other. Although this pore is not wide enough to accommodate the siderophore, and its accessibility is restricted by close contacts between S2, S5, and LC, the fact that it opens up toward the periplasmic faces has led to the proposal that it is a putative channel-forming region (Ferguson et al., 1998). However, any such opening is very much reduced in the ligand-loaded structure, as can be clearly seen in the figure, due to the unwinding of the switch-helix and the subsequent shifting of the N-terminal segment S1; in particular, residue Trp22, which in the siderophore-free conformation packs against S3, becomes buried in the entrance of the putative pore, reducing significantly its dimensions.
Water permeation through FhuA
To further investigate the existence of pore-forming regions in FhuA, we have analyzed the extent to which water permeates the barrel interior during the simulations, in both the ligand-free and bound states. Accessing this aspect of protein functionality is one of the strengths of the molecular dynamics simulation technique, as has been illustrated in studies on integral membrane proteins such as the human aquaporin AQP1 (de Groot and Grubmüller, 2001; Tajkhorshid et al., 2002), the glycerol transporter GlpF (Jensen et al., 2001), and the OmpF porin (Tieleman and Berendsen, 1998), as well as in the case of channel-forming peptides such as gramicidin (Chiu et al., 1999a,b; de Groot et al., 2002).
We have evaluated the distributions of the maximum residence time of water molecules in slices across the barrel domain (Fig. 9 A). The profiles are qualitatively similar in both simulations; in both the periplasmic space (z < 25 Å) and the extracellular region (z > 75 Å), most of the resident waters (around 80%) leave the corresponding slab in <1% of the 8 ns period analyzed. In contrast, in the central region (z 45 Å), a significant proportion of the resident waters remain in the region for 99% or more of the same period. The total number of water molecules resident in this slab in the period analyzed (from 2–10 ns) was 172 for FHUA and 125 for FHUA+F, respectively, reflecting the lower accessibility of this region in the siderophore-bound state. However, it is worth noting that in FHUA+F, the percentage of long-residing water molecules is more that two-fold than in the FHUA simulation (maximum values of 21.6 and 18.6% of the water molecules in the central slices compared to 8.4 and 7.6%, respectively). The existence of water molecules residing permanently in the plug-barrel interface is further illustrated by plotting the accessibility of each slice (defined as the proportion of the waters in the system that are resident in a given slice) as well as the water diffusion coefficient along the barrel axis (Fig. 9, B and C). As can be seen, the region 35 
z 50 Å is visited in FHUA by a small proportion of the water molecules in the system (12.3%), and by an even smaller proportion in FHUA+F (6.7%). These data correlate with a marked reduction in the diffusion coefficient with respect to bulk, again occurring to a greater extent in the FHUA+F simulation.
|
50 Å slab (Fig. 10). Trajectories are shown for five different intervals of residence times, illustrating the distinctive behavior of water at the barrel-plug interface. In particular, it is worth noting that several permeation events can be seen in the FHUA simulation. Visual inspection of the pathway followed by these water molecules reveals that they transverse the protein along the putative channel lined by regions S2, S5, and LC. However, permeation is abolished in FHUA+F in the timescale of the simulation. Here this channel is also solvated, although to a lesser extent due to the presence of residues of segment S1 at the entrance of the cavity.
|
; Karplus and McCammon, 2002). In all such studies, adequate sampling of the conformational space accessible for the species of interest is fundamental. For instance, in the context of protein dynamics, an optimal conformational sampling is most desirable if aspects of the protein functionality are to be understood on the basis of the motions observed during a simulation. In this section, we have addressed this issue by analyzing the dependence of the simulation B-factors on the length of the sampling window (see Methods for further details). As illustrated in Fig. 11 A, this simple assessment shows that the extent of the atomic fluctuations, averaged over the ß-barrel and plug domains, increases if longer windows of simulation are considered, and consequently, that the characteristic dynamic properties of the FhuA protein (as modeled in the simulation) have not been fully sampled. Specifically, this lack of convergence is more pronounced (that is, larger gradient in the plot) for the solvent-exposed domain of the protein than for the transmembrane region (Fig. 11 B), and is especially marked in the case of loop L8 (in contrast to, e.g., the S3 and S4 segments, Fig. 11 C).
|
; Im and Roux, 2002). However, it is often argued that the offset between crystallographic and simulation values is mainly due to the crystal static disorder. Although this is doubtless an important factor, our analysis indicates that insufficient conformational sampling in the simulations may account for a substantial part of this difference. |
DISCUSSION AND CONCLUSIONS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSION AND CONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
Comparison of the conformations adopted by the plug domain during the simulations, and in particular by segments S3 and S4, indicates that the subtle, possibly functionally relevant differences between the ligand-free and loaded states that were present in the crystal structures are preserved upon the introduction of thermal fluctuations. In the ß-barrel domain, the most noticeable conformational changes occur in the highly flexible, extracellular loop L8, whose motion seems to be dependent on the binding state of the protein. Specifically, it appears that in the presence of the bound siderophore, L8 changes its conformation relative to the transmembrane domain, so as to close the binding site, in a similar manner as in the structurally homologous ferric-citrate receptor FecA.
The other functionally important result from these simulations concerns the dynamics of water molecules within FhuA. The presence of the plug in the interior of the ß-barrel domain, and consequently, the numerous polar and nonpolar contacts between these domains, dramatically reduces the permeability of FhuA by water. This impermeability is specially marked in ferrichrome-bound FhuA, where the accessibility of water to the interior of the barrel domain is further reduced as a consequence of the unfolding of the switch helix and the subsequent displacement of the N-terminal segment. These findings support the view that a significant conformational change in the plug domain of FhuA is required for the siderophore to either passively diffuse or be translocated into the periplasm.
To extend the biological relevance of these and related simulations, further work is required to elucidate the possible path of conformational changes between the FHUA and the FHUA+F structures. It would be interesting to investigate possible conformational changes of the plug domain, including the feasibility of the hypothesis that it can be removed from the barrel. However, such a study might require an order-of-magnitude increase in computational power, even if an external force (i.e., steered MD (Lu and Schulten, 1999; Bockmann and Grubmüller, 2002)) was used to drive the conformational change.
How do these simulation results compare with structural biology and related experimental data? The recently determined structure of FecA (Ferguson et al., 2002), the TonB-dependent receptor of ferric citrate in E. coli, reveals that that two of the extracellular loops of the barrel, L7 and L8, change their conformation upon binding of the siderophore, seemingly so as to restrict the solvent accessibility of the binding pocket; this agrees with the simulation results for L8 in FhuA. Overall, our studies support the view that TonB-dependent receptors function in a similar way to air locks (Buchanan, 1999): ligand recognition and binding lead to the closure of a first hatch formed by the extracellular loops; subsequent structural changes in the periplasmic side of the protein, presumably associated with the unfolding of the switch-helix, might enhance interactions of the receptor with the Ton complex, which in turn would induce further conformational changes in the plug, allowing the ligand to enter the second hatch at the periplasmic space of the protein.
Assessment of the results
Molecular dynamics simulations have contributed significantly to the understanding of protein structure-function relationships (Hansson et al., 2002) as well as to the characterization of ligand-binding processes (Kollman, 1994; Åqvist et al., 2002). Similarly, MD has enabled an atomic-level description the interactions of channel-forming proteins with the solutes whose transport they mediate (Guidoni et al., 2000; Shrivastava and Sansom, 2000; Bernèche and Roux, 2001; de Groot and Grubmüller, 2001; Jensen et al., 2001; Tajkhorshid et al., 2002). In the context of membrane protein dynamics, simulations of the OmpF trimer, another outer membrane ß-barrel protein (Tieleman and Berendsen, 1998; Im and Roux, 2002) have provided valuable information on the role of loop L3 in conferring selectivity to this porin as well as on the ion permeation pathways. Simulations of the mechanosensitive channel MscL and of the bacterial OmpA have been used to study the relative flexibility throughout the protein, as a means by which insights into the mechanism of gating can be gained (Elmore and Dougherty, 2001; Gullingsgrud et al., 2001; Bond et al., 2002). On this basis, it seems that MD simulations should be able to contribute to the understanding of the structural, dynamical, and functional properties of TonB-dependent receptors.
However, it is important to consider the methodological limitations of current atomistic simulations of proteins. In addition to the simplifications inherent to most biomolecular force fields (notably, lack of electronic polarizability), it is still a matter of concern how to best represent the electrostatic interactions as well as the thermodynamical ensemble adopted. Similarly, we have discussed elsewhere (Faraldo-Gómez et al., 2002) the question of what might be the optimal way in which to embed a complex membrane protein within a lipid bilayer. In fact, the effect of the lipid composition of the bilayer model on the dynamics of membrane proteins remains to be determined. In the current work, a preequilibrated DMPC bilayer from a previous study (unpublished data) was used to reduce computational time. However, it is worth noting that this bilayer model differs significantly from the lipopolysaccharide-phospholipid asymmetric membrane in which FhuA naturally resides, and thus further work will be needed to draw general conclusions.
There remains the issue of how best to add water molecules to such a complex transporter, above and beyond those present with the x-ray structure (Yancopoulos et al., 2001; Yancopoulos and Mezei, 2002). On the basis of our simulation results, we speculate that reduced accessibility to the putative channel might well account for the absence of water permeation events. However, it should be noted that the extent to which the initial degree of solvation of this region might affect the permeation of water has yet to be explored. Similarly, it would be desirable to investigate the influence that the arrangement of H-bonds and salt bridges has on the water permeation process, e.g., by varying the protonation states of charged residues along the channel. This effect has been illustrated, for example, in recent simulations of the bacterial outer membrane protein OmpA (Bond et al., 2002).
Specifically, two aspects of the structural dynamics of FhuA have been the focus of the current study: first, the preservation of structural differences between the ligand-free and bound crystal forms of FhuA, specifically in S3 and S4 in the plug domain. Here, we have relied not only on the fact that the RMSD of the these segments with respect to the initial conformation is low and stable, but also on the persistence of the difference between their conformations relative to the ß-barrel, as indicated by their pairwise RMSD. This, combined with their relatively low dependence with the length of the sampling window (Fig. 11 B), particularly in S3, suggests that these truly represent two distinct states of the plug domain.
A second focus of our studies has been the motion of the barrel loop L8 and its possible role in closing the binding site and/or signaling the presence of the siderophore across the membrane. In contrast to S3 and S4, the conformational space of this slowly moving loop is poorly sampled. This is expected, given that its structure drifts consistently during the simulations away from the initial conformation. We note that the direction of the motion appears to be dependent on the presence or absence of the ligand. However, it cannot be fully excluded that, on longer timescales, a common conformational space may be shared by this loop in both binding states.
Finally, we support the view that incomplete sampling has to be added to the reasons for the generalized discrepancy between simulation and crystallographic B-factors, as has been discussed, for example, in Hüneberger et al. (1995). On the experimental side, errors may arise from, e.g., static disorder, as is apparent in the case of several membrane proteins; for instance, a comparison of the experimental B-factors for the topologically homologous TonB-dependent proteins FhuA (Ferguson et al., 1998) and FecA (Ferguson et al., 2002) yields average B-factors of 65 Å2 and 24 Å2, respectively. In view of these limitations on the parts of both simulation and experiment, it is desirable to improve the respective methodologies for assessing the dynamic properties of proteins.
|
ACKNOWLEDGEMENTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSION AND CONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
This work was funded by grants from the Engineering and Physical Sciences Research Council, The Wellcome Trust, The British Council, and La Caixa Foundation. The Oxford Supercomputing Centre provided computing facilities.
|
FOOTNOTES |
|---|
Submitted on September 12, 2002; accepted for publication April 15, 2003.
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSION AND CONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
Åqvist, J., V. B. Luzhkov, and O. Brandsdal. 2002. Ligand binding affinities from MD simulations. Acc. Chem. Res. 35:358–365.[Medline]
Bashford, D., and M. Karplus. 1990. pKa's of ionizable groups in proteins: atomic detail from a continuum electrostatic model. Biochemistry. 29:10219–10225.[Medline]
Berendsen, H. J. C., J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak. 1984. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:3684–3690.
Berendsen, H. J. C., D. van der Spoel, and R. van Drunen. 1995. GROMACS: a message-passing parallel molecular dynamics implementation. Comp. Phys. Comm. 95:43–56.
Bernèche, S., and B. Roux. 2000. Molecular dynamics of the KcsA K+ channel in a bilayer membrane. Biophys. J. 78:2900–2917.
Bernèche, S., and B. Roux. 2001. Energetics of ion conduction through the K+ channel. Nature. 414:73–76.[Medline]
Bockmann, R. A., and H. Grubmüller. 2002. Nanosecond molecular dynamics simulation of primary mechanical energy transfer steps in F1-ATP synthase. Nat. Struct. Bi
