INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Membrane proteins comprise ~30% of open reading frames (Wallin and von Heijne, 1998), and yet we know the x-ray structures of only ~40 such proteins (see, e.g., http://blanco.biomol.uci.edu/Membrane Proteins xtal.html). It is therefore essential that we extract the maximum possible information from the available structures. In particular, we wish to explore the dynamic behavior of membrane proteins, starting from the essentially static (time and space averaged) x-ray structures, as this is a key step in relating structure to biological function. It would also be valuable to be able to make predictions concerning membrane protein/bilayer interactions (Fyfe et al., 2001) on the basis of a crystal structure, often obtained in the absence of lipid molecules.

OmpA is a relatively simple bacterial outer membrane protein, for which two x-ray structures are available, at resolutions of 2.5 Å (Pautsch and Schulz, 1998) and 1.65 Å (Pautsch and Schulz, 2000). It is a member of the -barrel family of bacterial outer membrane proteins (Koebnik et al., 2000; Tamm et al., 2001). The N-terminal 171 amino acids (OmpA¹⁷¹) form a membrane-spanning domain comprising an eight-stranded anti-parallel -barrel. The C-terminal domain (the structure of which is not known) is thought to span the periplasmic space between the outer and inner membranes (Demot and Vanderleyden, 1994; Koebnik, 1995). The N-terminal domain has been described as forming an inverted micelle (Pautsch and Schulz, 1998) in that the exterior surface of the barrel is predominantly hydrophobic (as befits a membrane-spanning domain) whereas the interior of the barrel is formed by polar side chains and contains ~20 pore waters visible within the x-ray structure. These water molecules form three distinct groups; i.e., there is not a continuous water-filled pore running through the center of the crystal structure of OmpA¹⁷¹.

This latter observation is curious given several reports of formation of ion-permeable pores in lipid bilayers by OmpA¹⁷¹ (Sugawara and Nikaido, 1992, 1994; Saint et al., 1993). A recent study (Arora et al., 2000) has clearly demonstrated channel formation by the same N-terminal 171-residue fragment of OmpA that was used in the crystallographic studies. Furthermore, channel formation has also been seen for an orthologous protein OprF (Saint et al., 2000; Brinkman et al., 2000). Thus the question arises of how the x-ray structure of OmpA¹⁷¹ is related to the channel-forming properties of this protein. One hypothesis, which we will explore in this paper, is that a conformational transition must occur to convert the closed form of the pore seen in the x-ray structure to an open form that can form a continuous pore across a membrane. Some indication of the structural complexities of OmpA¹⁷¹ is provided by the recent NMR-derived structures of the N-terminal domain of this protein in dodecylphosphocholine (Arora et al., 2001) or dihexanoylphosphatidylcholine micelles (Fernandez et al., 2001). The studies of Arora et al. indicate that there is a gradient of flexibility along the axis of the -barrel on the micro- to millisecond timescale, which may relate to the fast, flickering conductance states observed when OmpA¹⁷¹ is reconstituted in planar lipid bilayers (Arora et al., 2000).

One way in which membrane protein structure dynamics may be explored is via molecular dynamics (MD) simulations. MD simulations have been widely employed to examine the conformational dynamics of pure lipid bilayers (Tieleman et al., 1997; Berger et al., 1997; Feller and Pastor, 1999; Schneider and Feller, 2001) and more recently have been used to investigate both peptides and proteins that span lipid bilayers (Forrest and Sansom, 2000). For example, MD simulations have been run for the following membrane proteins embedded in either a lipid bilayer or a bilayer-mimetic alkane slab: the porin OmpF (Tieleman and Berendsen, 1998), the ion channels KcsA (Guidoni et al., 2000; Shrivastava and Sansom, 2000; Bernèche and Roux, 2000, 2001) and MscL (Gullingsrud et al., 2001; Elmore and Dougherty, 2001), and the water-permeable pore aquaporin (and its homolog GlpF) (de Groot and Grubmüller, 2001; Zhu et al., 2001; Jensen et al., 2001). These simulations have revealed aspects of water and ion behavior within transbilayer pores (Sansom et al., 2000; Roux et al., 2000) and have provided information on protein/lipid interactions (Tieleman et al., 1999; Petrache et al., 2000) that complement those available from, e.g., crystallographic studies (Fyfe et al., 2001).

Here we report on a number of nanosecond simulations of the OmpA N-terminal domain embedded in either a dimyristoylphosphatidylcholine (DMPC) bilayer or in membrane-mimetic octane slabs. We explore the dynamics of protein and of water in these simulations. A combination of molecular modeling and MD simulation is used to propose a gating mechanism for the OmpA channel.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Protein model

Two OmpA¹⁷¹ crystal structures have been solved. In the 1.65-Å model (Protein Data Bank (PDB) code 1qjp) (Pautsch and Schulz, 2000), 34 loop residues could not be assigned a defined conformation. In comparison, the 2.5-Å model (PDB code 1bxw) (Pautsch and Schulz, 1998) has completely-described atomic positions for all loop residues, despite poor electron density. The C root mean squared deviation (RMSD) of the high- and low-resolution barrel regions was <0.1 Å, confirming that the two structures are very similar and that 2.5-Å OmpA¹⁷¹ is an adequate starting model. Before any simulations were performed, pK_A calculations were performed using the program UHBD combined with locally written code. The program UHBD (Davis et al., 1991), which solves the linearized Poisson-Boltzmann equation using a finite difference method, was used to calculate the difference in energy between each isolated amino acid and the corresponding residue in the protein-membrane environment, in both its protonated and unprotonated state (Bashford and Karplus, 1990; Adcock et al., 1998). The result was a series of titration curves for each ionizable residue, which allowed assignment of charges in the initial model. The resulting model was neutral overall.

Alternative solvation states of the initial model were investigated. Thirty-nine crystal waters were localized in the asymmetric unit. Twenty-seven of these were approximately categorized as lying within the bounds of the protein surface and were retained in the starting model for all but one of the following MD simulations. In one simulation, the program MMC (Mezei et al., 1985) (see http://fulcrum.physbio.mssm.edu/~mezei/mmc/) was used to perform a Monte Carlo (MC)-based solvation. Here, OmpA¹⁷¹ and the 27 crystal waters were embedded in a bilayer-mimetic slab of dummy atoms, and the system was then solvated with a 46 Å × 43 Å × 69 Å box of pre-equilibrated SPC water molecules. A Monte Carlo simulation in the grand canonical ensemble was then carried out, using GROMOS parameters at a temperature of 310 K. The protein and slab were treated as rigid, fixed solutes. Solvent insertion was cavity-biased at 100 × 100 × 100 grid points and force-biased with a lambda factor of 0.5. Each solvent shift was described by a 0.3-Å translation and a 30° rotation. After a million steps, the resultant solvent coordinates were used to select a new set of pore waters, which were used as the basis for setting up a new starting model.

Simulation system setup

A number of simulations were performed, as summarized in Table 1. For the DMPC simulation, the initial OmpA¹⁷¹ model was embedded in a pre-equilibrated bilayer as described in Faraldo-Gómez et al. (2002). Briefly, the GRASP solvent-accessible molecular surface (Nicholls et al., 1993) of OmpA¹⁷¹ was used as a template to remove lipids and to perform a short steered MD simulation of a solvated DMPC membrane to generate a cavity into which the initial OmpA¹⁷¹ model could be inserted. The system underwent 2.5 ns of MD (during the first 0.5 ns of which the protein coordinates were restrained) to allow equilibration of bilayer/protein interactions. This resulted in a stable area-per-lipid headgroup value of 65-70 Ų. A 5-ns production MD simulation was then carried out.

For the octane simulations, the appropriate initial protein (plus crystallographic and/or MC-inserted water) model was inserted in an octane slab. Each system was then solvated using a simple procedure of superimposition of a pre-equilibrated box of SPC waters followed by removal of any waters too close to either protein, octane, or crystallographic and/or MC-inserted water molecules. Protein-restrained MD simulations were run for 150-200 ps to allow equilibration, and subsequently 2-ns production MD simulations were carried out.

Simulation details

All simulations presented in this report were conducted using the GROMACS v2.0 (Berendsen et al., 1995) MD simulation package (www.gromacs.org). An extended united atom version of the GROMOS96 force field was used (Hermans et al., 1984). Every system was energy minimized before MD, using <1000 steps of the steepest descent method, to relax any steric conflicts generated during setup. Each neutral system was solvated with SPC waters and sodium and chloride ions corresponding to 1 M NaCl. During restrained runs, the protein was harmonically restrained with a force constant of 1000 kJ mol¹ nm². Electrostatics were calculated using particle mesh Ewald (Darden et al., 1993) with a 9-Å cutoff for the real space calculation. A cutoff of 10 Å was used for van der Waals interactions. All simulations were performed at constant temperature, pressure, and number of particles. The temperatures of the protein, octane/DMPC, and solvent (including both crystal and bulk water molecules, along with ions) were each coupled separately, using the Berendsen thermostat (Berendsen et al., 1984), at 300 K for octane simulations, or 310 K for the DMPC simulation, with coupling constant _T = 0.1 ps. The pressure was coupled using the Berendsen algorithm at 1 bar with coupling constant _P = 1 ps. For the DMPC simulation, the compressibility was set to 4.5 × 10⁵ bar¹ in all box dimensions. For the octane simulations, the compressibility in the x and y directions was set to zero to maintain the geometry of the bilayer-mimetic slab. The timestep for integration was 2 fs, and coordinates and velocities were saved every 5 ps. The LINCS algorithm was used to restrain bond lengths (Hess et al., 1997). The open state of OmpA¹⁷¹ was modeled using Quanta (Accelrys, San Diego, CA), and distance restraints were applied between the ionizable atom pairs of Arg138 and Glu128, with _disre = 10 ps with a force constant of 1000 kJ mol¹ nm².

Simulations were performed on a Linux workstation, an eight-node Beowulf cluster, or an SGI Origin 2000 using either four or eight parallel 195-MHz R10000 processors. Pore radius profiles and derived conductance estimations were calculated using HOLE (Smart et al., 1993, 1996, 1997, 1998). Secondary structure analysis used DSSP (Kabsch and Sander, 1983). Other analyses used GROMACS and/or locally written code. Molecular graphics was carried out using Molscript (Kraulis, 1991) and Povray (http://www.povray.org).

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Simulations

Eight different simulation systems for OmpA¹⁷¹ were set up (Table 1), seven with the protein embedded in an octane slab and one with OmpA¹⁷¹ in a lipid bilayer (DMPC). In all simulations, OmpA¹⁷¹ was placed such that the center of the membrane-spanning region, as defined by the midpoint of the two bands of Trp and Tyr side chains on the surface of the barrel (Fig. 1 A), was coincident with the center of the slab (Fig. 1 B) or bilayer (Fig. 1 C). In an attempt to obtain the most stable system, and to investigate the effects of different hydrophobic environments on protein dynamics, simulations with a DMPC bilayer and with differing octane slab thickness were compared. The mean distance between the two aromatic bands is ~21 Å. Hence, octane slabs of 20, 25, and 30 Å thickness were employed (Small_Oct, Med_Oct, and Big_Oct). The long axis of the barrel was oriented parallel to the z axis, where z is the normal to the octane slab or bilayer. Additionally, for the 25-Å slab three different orientations of the OmpA¹⁷¹ molecule were tested. Thus the OmpA¹⁷¹ molecule in Med_Oct (Table 1) was rotated by 90° about the z axis to give Med_Oct_90, whereas the OmpA¹⁷¹ molecule was tilted by 30° about a perpendicular to the z axis to give Med_Oct_Tilt. As will be discussed below, the aim of these three simulations was to explore the sensitivity of the results to the initial orientation of OmpA¹⁷¹.

View larger version (55K): [in this window] [in a new window]   FIGURE 1   (A) Schematic diagram of the OmpA N-terminal domain (drawn using VMD (Humphrey et al., 1996, and Povray (http://www.povray.org)) showing the two bands of amphipathic aromatic side chains (i.e., Trp and Tyr) and the presumed location of the membrane (horizontal gray lines). The extracellular loops (L1 to L4) are labeled. (B and C) Snapshots of simulation systems: Med_Oct (B) and DMPC (C) (drawn using Molscript (Kraulis, 1991) and Povray), showing the protein (in ribbons format colored from blue (N-terminus) to red (C-terminus)) and the octane slab (black) or DMPC bilayer (fatty acid tails in black; PC headgroups in red/blue/pink). Water is not shown for clarity. The L4 loop is labeled.

We also attempted to explore the influence of two different approaches to initial solvation of the OmpA¹⁷¹ barrel. In all but one simulation the crystallographic waters within the barrel were retained. In the simulation MC the crystallographic waters were supplemented by using a Monte Carlo search procedure (see above).

Structural drift and fluctuations

It is of interest to compare the drift of the OmpA¹⁷¹ structure from the initial, x-ray structure for the various simulation conditions. This can be assessed by calculating the RMSD for C atoms as a function of time. For most of the simulations in an octane slab, the all-residue C RMSD rose steadily before reaching a plateau of ~4 Å after ~0.5 ns (data not shown). This is in contrast to the plateau value of ~1.5 Å reached after 2 ns of the 5-ns DMPC simulation (Fig. 2 A), which in turn is more consistent with an ~2-Å C RMSD found in simulations of OmpF in a palmitoyloleoylphosphatidylcholine bilayer (Tieleman and Berendsen, 1998). Only Med_Oct and Med_Oct_90 gave comparable plateau values, each with 1.5-2-Å C RMSD after ~0.5 ns. Together, this suggests that the bilayer-mimetic octane slab allows more structural drift, on a 2-ns timescale, than does a phospholipid bilayer. However, with carefully tailored conditions, it seems that the membrane-mimetic system is capable of producing an equilibrated protein conformation of comparable divergence from the x-ray structure to that given in a phospholipid bilayer-based simulation, but on a considerably shorter timescale.

View larger version (34K): [in this window] [in a new window]   FIGURE 2   C RMSDs as a function of time. (A) All residues (solid black line) and -barrel residues (solid gray line) of the DMPC simulation; (B) Barrel residues of Med_Oct (solid black line), Big_Oct (solid gray line), and Small_Oct (dotted black line).

The large extracellular loops (L1 to L4; see Fig. 1 A) were expected to be more mobile; moreover, electron density for parts of these loops was not observed in either x-ray structure. Therefore we also analyzed C RMSD values for the -barrel residues only. These showed plateau values of <1.5 Å in all the simulations (Table 1), consistent with there being little loss of -sheet secondary structure (data not shown). Differences in -barrel C RMSD variation between simulations seem likely to be due to differences in the simulation environment. The Small_Oct and Big_Oct simulations produced the highest RMSDs of 1.2-1.4 Å (Fig. 2 B). In comparison, the barrel RMSDs for the DMPC (Fig. 2 A) and Med_Oct simulations (Table 1 and, e.g., Fig. 2 B) were much lower at 0.7-0.9 Å. This suggests that, on a nanosecond timescale, simulations using the medium octane slab yield comparable stability of the transmembrane domain to the more realistic (and less fluid) model of the bilayer. Comparison with results from the NMR solution studies of Arora et al. (2001) is instructive. An alignment of the lowest-energy NMR structure (PDB code 1g90) and the 1.65-Å x-ray structure (Pautsch and Schulz, 2000) for the 108 best-fitting barrel and turn residues yielded an RMSD of 1.3 Å (Arora et al., 2001). An RMSD of 1.7 Å was achieved between the starting (2.5-Å) x-ray structure (Pautsch and Schulz, 1998) and the 5-ns simulated structure using the same 108 residues.

Analysis of the protein's secondary structure showed a high degree of barrel stability in all simulations (data not shown), consistent with the RMSD analysis. Even for Small_Oct, which displayed the greatest structural drift, very little loss of secondary structure occurred. Indeed, ~10 L4 residues converted to -sheet structure after ~0.8 ns, extending the two shortest -strands (7 and 8) of the crystal structure. Aside from this, in various simulations, the four longest -strands (3 to 6) were occasionally seen to transiently increase or decrease in length by 2-4 residues. These results agree qualitatively with the OmpA solution NMR studies (Arora et al., 2001): for the lowest-energy NMR structure, strands 7 and 8 are one residue shorter, and strands 3 to 6 are three residues shorter than those of the crystal structure.

Local protein mobility was analyzed by calculating the time-averaged RMS fluctuation of each residue (i.e., C RMSF) and converting them to the equivalents of crystallographic B-values (Fig. 3); this was performed for the equilibrated time period of each simulation (the final 1.5 ns or 3 ns for all octane simulations or for the DMPC simulation, respectively) to remove coupling between total RMSDs and calculated atom fluctuations. Only -barrel C residues were used during the fitting process to the initial, low-resolution x-ray structure, to ensure that the low-mobility barrel residues did not artifactually attain higher motion due to the loop residues.

View larger version (33K): [in this window] [in a new window]   FIGURE 3   Calculated B-values (temperature factors) for C atoms plotted as a function of residue number. (A) The B-values from the experimental high-resolution crystallographic structure (thin broken black line) (Pautsch and Schulz, 2000) are compared with those derived from the RMS fluctuations of the C atoms during the DMPC simulation (thick solid gray line). (B) The calculated B-values for the Small_Oct (thin solid black line), Med_Oct (thick solid gray line) and Big_Oct (thin broken black line) are compared. The labels, indicating the regions with high B-values, are for the extracellular loops (L) and periplasmic turns (T).

In all simulations, there was a clear correlation between mobility and local structure, with small, medium, and large B-values for the barrel strands, the periplasmic turns, and the extracellular loops, respectively. Despite performing the RMSF calculations using only the simulation trajectory periods judged to have become equilibrated (according to RMSD analysis), there was still a correlation between the plateau RMSD values and overall atom mobility. Thus, of all simulations, Small_Oct and Big_Oct generally had the highest temperature factors. In comparison, Med_Oct and DMPC both demonstrated lower atom fluctuations on average. Possible reasons for this difference in mobility are discussed in the following section.

The calculated simulation temperature factors were compared with the 1.65-Å crystal structure (Pautsch and Schulz, 2000) values to judge the extent of correlation of residue-by-residue mobility between the various simulations and the x-ray structure. Of course, it is not possible to compare B-values for the high-mobility loops, for which there is no crystallographic estimate. For the rest of the protein, despite the qualitative agreement according to secondary structural features, simulation B-values are consistently lower than those of the crystal structure. For comparison, a number of studies of soluble proteins have revealed discrepancies between B-factors calculated from MD simulations and those from crystallography. Studies on bovine pancreatic trypsin inhibitor and lysozyme (Hunenberger et al., 1995), and also on crambin (Caves et al., 1998), suggest as one explanation for such differences to be incomplete conformational sampling during individual simulations. Here, however, none of the OmpA barrel C RMSDs between each possible pair of final (equilibrated) simulation structures significantly exceeds the largest trajectory plateau RMSD value of the respective pair. This suggests that a limited conformational landscape is available to the OmpA barrel under a variety of conditions (within a membrane context), consistent with the high stability of -barrels observed experimentally.

Additionally, differences between the crystal and membrane environments will probably lead to differences in mobility. It is hard to estimate the extent of this effect because there are probably more detergent molecules present in the crystal than are observed, nestling in inter-protein grooves within the endless OmpA¹⁷¹ fiber formed along the crystal c-axis of the high-resolution structure (Pautsch and Schulz, 2000). However, intuitively, one might envisage that a lipid membrane environment would more tightly constrain the OmpA barrel. Perhaps only an MD simulation of OmpA in its crystalline lattice could answer this question. Finally, one should consider limitations in techniques. In the case of crystallography, some of the apparent protein mobility may be contributed by (static) crystal disorder. In the case of MD, the force field may lead to unrealistic dynamics. Unfortunately, these factors are difficult to assess.

For comparison, the C RMSDs for each residue in the 10 lowest-energy NMR structures (Arora et al., 2001) were calculated and converted to B-values, using the average coordinates of these as a reference for fitting. In all cases, the B-values estimated from the NMR solution studies were significantly larger than those of the MD simulations presented here. However, it is impossible to know how much of this apparent increase in protein mobility results from the detergent micelle environment, as the estimated B-values are probably artificially high due to the lack of sufficient NOE data, particularly in the loop regions.

Protein-environment interactions

Simulations with different thickness octane slabs reveal how OmpA¹⁷¹ interacts with a bilayer-mimetic environment, for which the interfacial region (where octane and water molecules are intermixed) is ~5 Å thick. Despite the presence of some degree of protein/octane hydrophobic mismatch there is no evident distortion of the octane slab to match the spacing (~21 Å) of the two aromatic bands on the surface of OmpA¹⁷¹. In all cases, the lower (i.e., periplasmic) aromatic belt seems to be more closely localized within the interfacial region compared with the upper (i.e., extracellular) belt. This may correlate with a somewhat larger and more regular amphipathic aromatic belt (i.e., surface-exposed Trp and Tyr side chains) at the lower (periplasmic) end of the barrel (Fig. 1 A).

The Med_Oct slab size was judged to be optimal, because the peaks of aromatic belt density lie at the water-octane interfaces (Fig. 4 A). For the DMPC simulation (Fig. 4 B), although the maximum width of hydrophobic density is larger than that of Med_Oct, there is also a deeper penetration of polar atoms within this density, due to the presence of lipid phosphates here (and a subsequent reduction in the energy barrier required for water to penetrate). Therefore, the belt aromatic side chains still lie at the polar-hydrophobic interface, in the region where solvent, lipid polar headgroups and apolar tails all coincide. Once again, the lower aromatic belt is more closely associated with the interfacial region than the upper belt. The interfacial positioning of the aromatic belts in both the Med_Oct and DMPC simulations is consistent with the similarly low drifts and high stability of the systems, as discussed previously.

View larger version (15K): [in this window] [in a new window]   FIGURE 4   Trajectory-averaged system density as a function of the z coordinate for Med_Oct (A) (with arrows representing variation in width of octane slab density for Small_Oct and Big_Oct) and the DMPC simulation (B). The protein aromatic belts (black filled areas) can be seen to localize at the interface between the octane (A) (black dashed line)/water (solid gray line) interface, or the lipid tail (B) (black dashed line)/polar headgroup (gray dashed line) and water (solid gray line) interface. The lower aromatic belt (at z  25 Å) localizes better with the polar/hydrophobic interfaces in each case, whereas the overall protein density (solid black line) is preferentially shifted toward the extracellular side.

Counting the numbers of H-bonds between the aromatic side chains and polar groups (i.e., water molecules and lipid headgroups) reveals that, consistent with the preferential localization of the lower belt at the interface, the number of H-bonds of the lower belt is always substantially greater than that of the upper belt (Table 2). In terms of individual simulations, those using a medium-sized octane slab have the greatest number of H-bonds, presumably because the best balance is achieved between barrel surface area buried in octane and polar sections of aromatic belt side chains exposed to solvent. Perhaps counterintuitively, the aromatic belt side chains in the DMPC simulation make somewhat fewer H-bonds to polar groups in comparison with the octane simulations. However, there is also less fluctuation in the number of H-bonds on the picosecond timescale over the course of the simulation. It appears that more long-lived H-bonds between the protein and the phosphate headgroups (and structured waters) compensate the entropic loss, in comparison with the exclusively water-mediated H-bonds in the octane simulations.

As mentioned above, no distortion of the octane slab to match the spacing of the aromatic bands occurs. However, the OmpA¹⁷¹ barrel does tend to tilt relative to the bilayer/slab normal on an ~200-ps timescale, and this may be a method by which the protein optimizes aromatic band positioning and maximizes the buried hydrophobic surface area. In the medium-octane simulations in which the barrel begins in an untilted state, it seems to stabilize quickly at an angle of 5-10°, consistent with it being the most stable octane simulation (Fig. 5 A). In comparison, in Small_Oct, the barrel tilts to a stable 15-20° angle after ~0.5 ns; the longer -strands on the extracellular side appear to be leaning toward the octane slab, increasing the total amount of buried, hydrophobic surface area (Fig. 5 A). In Big_Oct, the barrel rapidly fluctuates between 0° and 20° over 2 ns; this may represent an entropic conflict between maximizing buried surface area and exposing the polar loops to solvent (Fig. 5 A).

View larger version (32K): [in this window] [in a new window]   FIGURE 5   Barrel tilt angle as a function of time for simulations: (A) Med_Oct (solid black line), Big_Oct (solid gray line), and Small_Oct (dashed line); (B) Med_Tilt (solid line) and Med_90 (dashed line); (C) DMPC.

Unsurprisingly, the barrel tilt behavior of Med_Oct_90 is very similar to that of Med_Oct (Fig. 5 B). In contrast, in Med_Oct_Tilt, where the barrel has an initial angle of 35°, OmpA¹⁷¹ stabilizes to an angle of 20° within just 100 ps, followed by a fluctuation between 15° and 25° during the remainder of the simulation (Fig. 5 B). The protein could be trapped in a local energy minimum here; the 2-ns period may not be long enough for OmpA¹⁷¹ to equilibrate, considering its initially somewhat extreme orientation. Finally, in DMPC, consistent with its low RMSD and the less fluid nature of the bilayer model compared with octane, the barrel tilt angle fluctuates between just 0° and 10° over the entire 5 ns (Fig. 5 C). Overall, it is difficult to establish whether the tilting of the barrel in each simulation has reached equilibrium; it is certainly conceivable that a small amount of barrel tilting (~5-10°) is a general property on the nanosecond timescale.

Recently, an empirical Monte Carlo minimization procedure, IMPALA (Ducarme et al., 1998), has been developed to predict the depth and angle of insertion of several membrane proteins using a simple empirical energy function to describe a lipid bilayer. When applied to OmpA¹⁷¹, an energy minimum was found corresponding to a shift of the protein mass center by 7 Å (in an extracellular direction) from the bilayer center, and an angle of insertion of reference strand 2 of 52° with respect to the membrane surface (Basyn et al., 2001); this angle is equivalent to a barrel axis tilt relative to the bilayer normal of ~9°. We have estimated equivalent average values for the stable time periods in each trajectory (i.e., the last 1.5 ns in all octane simulations and the last 3 ns in DMPC). In our most stable simulations, there is relatively good agreement with the IMPALA method; thus, Med_Oct gave a barrel displacement of 6.3 ± 0.6 Å and tilt of 9.0 ± 1.6°, and DMPC gave a displacement of 6.9 ± 0.6 Å and a tilt of 5.1 ± 2.4° (Fig. 5 C). Of course, the IMPALA approach treats the membrane and protein as completely rigid bodies and uses an empirical bilayer energy function, but nevertheless, our results suggest that there may be some benefit in applying such a procedure while setting up initial models for full atomistic MD simulations.

Pore properties

A Monte Carlo solvation procedure was used to investigate alternative solvation states of the initial model. Fairly good agreement was achieved for Monte Carlo solvent molecules within the membrane-embedded region and the crystal waters, suggesting that the observed crystal-form cavities were already fully solvated (Fig. 6). However, toward the outer regions of the barrel, extra waters were introduced. These may have important effects on global changes in protein structure and may induce expansion of the -barrel domain.

View larger version (25K): [in this window] [in a new window]   FIGURE 6   Histograms comparing the numbers of waters within the pore of OmpA171 for the original crystallographic waters (dashed line) and that produced using a Monte Carlo solvation procedure (solid line). A sampling width of 2.5 Å was used for the pore radius z coordinate (x axis). Good agreement is observed within the center of the pore, but additional water can be observed toward the pore mouths in the MC model.

Analysis of pore profiles during the simulation trajectories highlights the dynamic nature of the -barrel interior. Even for the most stable Med_Oct simulation, the internal volume fluctuates by up to ~100 ų on a 10-ps timescale, and over the course of the simulation, it varies between a minimum of 100 Å³ (~3 waters) and a maximum of 600 Å³ (~20 waters) (data not shown). There is always some degree of expansion of the pore, particularly toward the extracellular and periplasmic mouths of the -barrel. Thus, on the extracellular side of the barrel, all pore radii measurements in the region of the ionizable residues Glu140, Lys12, and Arg96, or polar ring, have standard deviations that are higher than equivalent crystal structure values. Pautsch and Schulz (1998) proposed that an N-terminal periplasmic cover formed by the residues 1-4, stabilized by a H-bond between Asp92 and the amide of Ala1, would block the periplasmic barrel entrance. Indeed, such a H-bond was maintained in each simulation; however, residues 1-4 appear to be relatively mobile, and the region repositions itself to create a wider periplasmic mouth. Therefore, in all cases, this leaves the Arg138-Glu52 gate as the major constriction of the potential pore (see Fig. 7 A).

View larger version (24K): [in this window] [in a new window]   FIGURE 7   Schematic diagram (produced using ISIS draw (MDL Information Systems, San Leandro, CA) showing the location of key constriction side chains within the pore of OmpA171. In the crystal structure (A) several ionizable side chains appear to prevent permeation by water or ions, in particular the Arg138-Glu52 salt bridge, or gate, in the center of the pore. In our model of the open state (B) a salt bridge switch has occurred such that Arg138 now is paired with Glu128. Additionally, the polar ring and periplasmic cover are proposed to be sufficiently mobile to allow ion and/or water permeation, as suggested by the arrows.

The pore radius profiles for the OmpA¹⁷¹ starting structure and the Med_Oct simulation are shown in Fig. 8, A and B; other octane simulations demonstrated similar internal fluctuations, although in Small_Oct and Big_Oct, OmpA¹⁷¹ showed more fluctuation at the barrel openings, consistent with the higher RMSDs. The DMPC simulation unsurprisingly showed lower pore radius mobility; the timescale sampled was almost certainly too short to observe such dynamic conformational shifts.

View larger version (22K): [in this window] [in a new window]   FIGURE 8   Pore radius profiles (produced using HOLE (Smart et al., 1996)). (A) The pore profile for the starting model of OmpA171 is shown, and the major constrictions (pc, periplasmic cover; g, gate; pr, polar ring) to permeation are indicated. (B) The trajectory-averaged pore profile with standard deviation bars (gray) for the most stable octane simulation, Med_Oct, reveals that only around the Arg138-Glu52 gate region is the pore radius unlikely to expand to the minimum dimensions required for a water molecule to pass. (C) The trajectory-average pore profile with standard deviation bars for the simulation of the modeled OmpA171 open state is shown. Notice that now the pore radius around the gate region is able to expand beyond 1.15 Å.

Modeling the open state

Regarding the evidently dynamic nature of OmpA¹⁷¹, one may ask whether a water molecule could pass through the entirety of the pore. If a water molecule is unable to pass through the pore, it is unlikely that an ion could. Even in the most stable Med_Oct simulation, only at one point along the pore axis is the radius unable to transiently reach a size large enough (>1.15 Å) to accommodate a water molecule (Fig. 8 B). This region is equivalent to the Arg138-Glu52 gate. Consistent with this, z (pore) axis trajectories of individual water molecules demonstrate that in the more mobile Small_Oct, Big_Oct, and Med_Oct_Tilt simulations, water molecules can be seen to cross-the polar ring. Moreover, in all simulations, water crossed the periplasmic cover (Fig. 9 A). In contrast, in none of the simulations do water molecules cross the gate.

View larger version (36K): [in this window] [in a new window]   FIGURE 9   Trajectories of selected water molecules projected onto the z (pore) axis (solid black lines) with respect to the protein aromatic belts (solid gray lines). Note that for clarity the majority of water trajectories are not shown. (A) The water trajectory plot for Big_Oct reveals that water is able to pass both the polar ring (pr) and periplasmic cover (pc), but no waters are able to pass the gate (g). (These major constriction regions are indicated by arrows.) (B) The water trajectory plot for the simulation of the modeled open state reveals that water molecules are now able to pass the Arg138-Glu52 constriction (g). One water molecule (broken line) is observed to traverse the entire length of the pore.

The dynamic state of the pore demonstrated here, combined with a wider body of experimental evidence, suggests the possibility that the hydrophilic interior of the protein may be able to adopt a subtly different conformation from that observed in the crystal structure. Analysis of the internal H-bonding network reveals the existence of the carboxylate group of Glu128 lying ~2 Å from one of the guanidinium NH₂ groups of Arg138. (Although Glu128 lies opposite Lys82, the carboxylate and amine groups are separated by over 4 Å; moreover, water is able to freely pass between these two residues during MD; Fig. 10 A). This suggests the possibility that the Arg138 side chain may be able to form alternative salt bridges with either Glu128 or Glu52 (Fig. 7, A and B). Indeed, during simulation, the side chain of Arg138 displayed significant conformational mobility, such that the guanidinium N proton switched from pointing toward the extracellular side of the pore (as in the crystal structure) to an orientation in which it pointed toward Glu128. To test this hypothesis, a model was constructed in which the Arg138 side chain was isomerized to a new rotamer (Dunbrack and Karplus, 1993), orientating it toward Glu128 (Fig. 10 B). For the normal and re-modeled systems, the sum of short-range Coulombic interaction energies around the gate region were approximately equal (data not shown). Thus the new conformation has about the same potential energy as that observed in the x-ray structure. A 2-ns simulation (open) was performed, during which a restraint was applied between these two side chains. The results revealed that such a conformational change enables water-crossing events to occur past the gate region, as revealed by comparison of water trajectories (Fig. 9 B). The pore profile for this MD trajectory confirms that an increase in radius around this region has indeed occurred (Fig. 8 C).

View larger version (54K): [in this window] [in a new window]   FIGURE 10   System snapshots of the gate region: (A) a standard simulation (Med_Oct), in which water is unable to pass the Arg138-Glu52 salt bridge (indicated by broken black lines); (B) modeled open simulation, in which Arg138 is now paired with Glu128 (broken black lines); note that a single file of water molecules is now present.

It is possible to predict the approximate conductance of a pore by treating it as a simple Ohmic resistance equivalent to a cylinder of electrolyte with reduced ionic mobility (Smart et al., 1997). Briefly, the approximation of the channel conductance assumes that a pore of given dimensions is filled with an electrolytic solution of defined resistivity. Additionally, because diffusion of ions is somewhat slower within a pore of molecular dimensions (Smith and Sansom, 1998, 1999) the conductivity of an ionic solution within a channel is smaller than that of bulk solution, an empirically-based correction factor dependent on the minimum radius of the channel is used (Smart et al., 1997). This method was used to predict the conductance of the open state of OmpA¹⁷¹ over the course of the simulation, assuming 1 M KCl as the electrolyte. As shown in Fig. 11 the predicted conductance fluctuates on an ~200-ps timescale between ~30 and 80 pS, giving an average conductance value of 51 ± 11 pS; values of ~ 60 pS (with a range from 50 to 80 pS depending on the exact protein construct used) are observed under similar experimental conditions (Arora et al., 2000). Thus, the model of the open state is predicted to have a conductance of the same magnitude as that seen experimentally for OmpA¹⁷¹.

View larger version (21K): [in this window] [in a new window]   FIGURE 11   HOLE-based (Smart et al., 1997, 1998) prediction of pore conductance as a function of time for the modeled open-state simulation, assuming 1 M KCl as the electrolyte. The average predicted conductance is 51 ± 11 pS, which is similar to the experimental value of ~60 pS (Arora et al., 2000) (both shown as dotted horizontal lines).

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Conclusions

The significant conclusions to be drawn are threefold. First, it is evident that starting configuration and membrane environment can influence the protein behavior during an MD simulation of OmpA and that it is possible to reproduce on a shorter timescale the characteristics of an explicit lipid bilayer model using a membrane mimetic. Second, our simulations demonstrate the dynamic behavior of a simple -barrel outer membrane protein, in the context of stable MD simulations. In all simulations, some degree of pore expansion occurred, particularly at the barrel mouths, consistent with recent NMR studies (Arora et al., 2001). Finally, a small perturbation of the Arg138-Glu52 gate enabled complete permeation by several water molecules and resulted in a calculated conductance very similar to experimental data. Thus, we hypothesize a gating mechanism of OmpA¹⁷¹, achieved through the rotamerization of Arg138, alternating between being paired to Glu52 and to Glu128 (Figs. 7, A and B).

It should be noted that the ~60-pS channel-opening events observed by single-channel conductance measurements for native OmpA occur on a millisecond timescale (Arora et al., 2000) whereas the physical process of our proposed gating mechanism may well occur on a much shorter timescale. Unfortunately, standard MD simulations do not allow us to address such a problem. As we have mentioned, the local electrostatic energy around the gate region is similar in both closed and open conformations, but one requires an estimate of the transition-state energy barrier between the two states to obtain an associated rate constant. Conceivably computational methods could be used to find a minimum energy path between the closed and open states (Smart, 1994; Smart and Goodfellow, 1995).

Nevertheless, we believe that the timescale of the driving force for the gating mechanism is the more relevant focus in relation to experimental measurements of channel opening and closing. The two forms of OmpA¹⁷¹ (closed and open) are proposed to be similar-energy conformations separated by a relatively high barrier. Also, the open state may be stabilized on a millisecond timescale by the presence of ions and/or a transmembrane potential. We note that a gating mechanism based on switching of ion pairs to open and close a channel is not unprecedented, as a similar mechanism has been proposed for the annexin family of amphipathic, calcium-dependent phospholipid-binding proteins, which form ion channels in vitro, in this case driven by voltage (Benz and Hofmann, 1997).

Beyond this, the main deficiencies of these simulations result from limitations in computing power; a total of 19 ns has been simulated here, yet lipid and protein motion is known to occur on a longer timescale. Statistically speaking, one would like to extend simulation times further, or run multiple, short simulations to sample as much conformational space as possible. This would ensure greater confidence that the observed water-crossing events and pore dynamics are a significant system property. Of course, explicit extrapolations to experimental data are not yet possible. Consider ion passage. The experimentally measured channel conductance g = 60 pS with V = 100 mV gives an electrical current I = gV = 6 × 10¹² A  3.7 × 10⁷ monovalent ions s¹. Thus, the mean time spent by a single ion in passing from one side of OmpA¹⁷¹ to the other is ~27 ns, a timescale not yet readily accessible to standard MD methods with system sizes of ~40,000 atoms. Similarly, the experimentally observed flickering current measurements (Arora et al., 2000), presumed to be equivalent to gating transitions between closed and open states of OmpA¹⁷¹, occur on a timescale of ~ 1 ms.

The latter timescale is the motivation for our attempt to model an open state and then to explore this state by nanosecond simulations. Our conductance prediction should be considered somewhat approximate, due to the short simulation time and the arbitrary element in the modeling procedure. Moreover, the conductance estimation procedure uses an empirical correction factor, so that the accuracy of such a prediction is limited (Smart et al., 1997, 1998). So how can one have any real certainty whether or not the perturbation of the Arg138-Glu52 salt bridge is equivalent to the real OmpA¹⁷¹ channel opening event? First, we have noted that it is the region past which water does not pass in any simulation. Obviously, if ions are to permeate, then so too must water. Second, our conductance estimations are extremely encouraging given that similar experimental values (~60 pS) have been achieved for reconstituted OmpA¹⁷¹ (Arora et al., 2000). Third, as mentioned already, a similar mechanism is proposed for the annexins. Finally, the -sheet-rich N-terminal domain of OprF has been shown to form channels of 360 pS in lipid bilayer membranes; a model based on the orthologous OmpA¹⁷¹ suggests that OprF may form larger conductance channels because of the lack of conservation of pore-constricting residues (Brinkman et al., 2000). Indeed, a pore profile of the model (data not shown) confirms that there is a larger radius around the region equivalent to the gate of OmpA¹⁷¹.

Future directions

The small size and high stability of OmpA, along with its amenability to study by various biophysical techniques, makes it an ideal subject for comprehensive MD studies and comparative analysis with experimental results. Beyond this, extension of the simulation study of the modeled open state may be useful. Free energy perturbation methods could be used to estimate the free energy differences between the Arg138-Glu52 and Arg138-Glu128 conformations, and nonequilibrium MD might be used to attempt to drag ions through the pore of OmpA¹⁷¹. Additionally, the membrane environments simulated here do not realistically represent the complex and heterogenous lipopolysaccharide (LPS) of the outer membrane, which may affect the local dynamics and electrostatics of embedded proteins. For example, phage activity via OmpA¹⁷¹ recognition has been shown to require LPS, and it has been proposed that LPS may be involved in modifications of pore formation by the ortholog OprF (Freulet-Marriere et al., 2000). Thus, attempts could be made to make simulations more realistic, by modeling OmpA¹⁷¹ in a LPS-containing environment (Lins and Straatsma, 2001). Finally, a potential experimental approach may be to mutate the residues we have hypothesized to act in the gating mechanism.