INTRODUCTION

TOP ABSTRACT INTRODUCTION THE MODEL SYSTEM RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The passage of ions through a membrane is a fundamental physical process that is common to all organisms. This process is the key mechanism in homeostasis, cell excitation, flagellar motion, ATP synthesis and, as recently indicated (Groves, 2000), may also cause local translational motion and demixing effects in biomembranes. The experimental systems for monitoring the ionic current are well established, and parameters like single channel conductance or ionic selectivity are readily measured. However, despite accumulated knowledge about channel proteins and their structure, our capacity to predict the conductance of channels on the basis of their structures is limited.

In principle, all the forces that operate on an ion, from the entry to the vestibule (Jordan, 1982) and during its passage, are derived from the fixed charges on the protein surface and their relative placement with respect to the dielectric boundary. These forces, as modulated by the protein dynamics and ionic atmosphere, should be sufficient to predict the passage rate of ions along a protein channel. On the basis of this information, the ionic conductivity of channels can be predicted.

Three methodologies that correlate the structure with the ion flux are currently presented: 1) the HOLE program (Smart et al., 1997), which is based on the geometry of the conducting channel. According to this procedure, the program calculates the maximum space within the pore, and the conductance is calculated according to the specific conductance of the bathing electrolyte. A correction factor had to be introduced to match the calculation with the observation by reducing the ionic mobility inside the channel. The lower mobility was attributed to the restricted motion of the solvent molecules near the protein's surface (Fischer et al., 2000). The prediction of the flux relies on the atomic-resolution structure but ignores the intra-channel electrostatic field and the dynamics of the protein. 2) The detailed Poisson-Nernst-Plank (Kurnikova et al., 1999). In this case the structure, charge distribution, and ionic atmosphere are all accounted for to generate the precise electrostatic potentials within the conducting channel using the Poisson-Boltzmann equation (Dennis et al., 2000; Sandberg and Edholm, 1997; Sitkoff et al., 1994; Polozov et al., 1999; Pellequer and Chen, 1997; Marino et al., 1995; Holst et al., 1994; Oberoi and Allewell, 1993). The steady-state ion flux is calculated by the application of the Debye-Smoluchowski equation. This method was successfully applied for reconstructing the conductance of the gramicidin channel, at the expense of heavy computation time (from 10 min to a few hours of computation time for each point in a grid of 10⁶ (Kurnikova et al., 1999)). At present, its application to more complex protein structures is not practical. 3) Brownian motion simulation of ions in a channel (Crozier et al., 2001; Jakobsson, 1998; Phale et al., 2001; Schirmer and Phale, 1999). This approach is based on the detailed electrostatic potential at the vestibule and inside the channel that was calculated by the Poisson-Boltzmann equation (the protein dynamics were not included in the calculations). The potential field accounted for the ionic screening through a continuum approximation. Within the defined space, molecular dynamics calculations were carried out, each time for a single particle (anion or cation) that was released at the vestibule and its random Brownian motion was followed up to 10⁷ steps. As many trajectories ended with the particle wandering back to the bulk, the calculations were repeated 5000 times for positive particles and an equal number of negative particles. Of the many calculated trajectories, only a small fraction (~10%) transversed the whole length of the channel. These "successful" trajectories were subjected to statistical analysis that functioned as a mathematical sieve that emphasized the repeated implementation of the Boltzmann exponential term (exp(E/kT)) that correlates the probability of motion in a given direction with the gradient of the free energy. This mode of calculation is very comprehensive, tracing both the ion in the channel and those wandering at the channel's vestibule. Because of that, the correlation between the probability of crossing the channel and the measured ionic selectivity and single channel conductance is high (Schirmer and Phale, 1999; Phale et al., 2001).

The simulation of Brownian motion trajectories is overburdened by heavy computational redundancy (Jakobsson, 1998), wasting most of the computation time on the unproductive attempts of ions to get out of traps or wander in "nonproductive directions." These futile calculations can be avoided by reversing the order of operations; first, to determine, on the basis of detailed electrostatic potential maps, what will be the most probable trajectory, and then to transport the particle along the selected trajectory. This algorithm (propagation along a structure-supported trajectory (PSST)) reduces the diffusion in a three-dimensional space to a one-dimensional diffusion problem with a computation time of a few seconds.

The present publication and the adjacent one implement the reversal of order in calculation of ion passage, using as a model the large pore channel family of ion conducting proteins (Benz et al., 1984, 1989; Berrier et al., 1997; Eppens et al., 1997b; Jap et al., 1991; Phale et al., 2001; Schirmer and Phale, 1999; Struyve et al., 1993; Sutcliffe et al., 1983; Watanabe et al., 1997). These proteins were selected because of the presence of ample information concerning both structure and single channel conductance (Eppens et al., 1997a; Jap et al., 1991; Koebnik et al., 2000; Li et al., 1998; Vangelder et al., 1996; Watanabe et al., 1997). Another advantage of the porin family is the large size of the pore, which ensures that the passage will neither be delayed by desolvation barriers nor affected by closure of the channel during structure fluctuations common to very narrow channels (Tieleman and Berendsen, 1998).

In the present article we describe an algorithm that selects, within the diffusion space, the most probable trajectory and we test its capacity to reconstruct time-resolved measurements of a single proton released inside the PhoE channel. The experimental system (Gutman et al., 1992b) consisted of a single pyranine molecule (8-hydroxypyrene-1,3,6-trisulfonate, OH) that was bound by electrostatic interactions inside a PhoE channel. A short laser pulse caused proton dissociation from the excited molecule, and the measured time-resolved fluorescence decay curve was reconstructed by the geminate recombination model of Agmon (Pines et al., 1988). In a previous study (Shimoni et al., 1993), the signals were analyzed using a simplified symmetric Coulombic potential field. In the present study the electrostatic potentials were calculated according to the Poisson-Boltzmann equation, where the input was the detailed structure of the protein, including the partial charges of all atoms. The intra-cavity electrostatic potential was mapped in thin slices perpendicular to the long axis of the channel, and the path that offered minimal resistance to the proton's propagation was selected. This procedure reduced the complexity of the system to a simple case of a linear array of sites, and the transition probability between adjacent sites is determined by the gradient of the electrostatic field and the appropriate entropy terms.

A search that was carried out within a limited parameter space yielded the terms that are sufficient to characterize the proton transfer reactions in the intra-cavity space: 1) the rate constants of the reversible proton transfer between the excited pyranine molecule and the adjacent water molecules; 2) the reactivity of the carboxylates lining the intra-cavity with free diffusing proton; and 3) the dielectric constant of the aqueous phase inside the channel (_int-cav  55). This value is the first experimental confirmation of Sansom's (Breed et al., 1996) estimation of the dielectric constant of water in -sheet barrel structure.

In the accompanying publication we tested the extent to which the algorithm, which had been refined to reconstruct the propagation of proton in the channel, can be used to predict the conductance of five porin proteins (PhoE, OmpF, and three general diffusion porins from R. blastica (WT and two mutants)). As will be shown, the electrostatic potential calculation based on  = 50 yielded values that predicted passage times compatible with those derived from the single channel conductance measured values of the three proteins.

Finally, the advantage of the present mode of calculation is discussed with respect to the computation time, which is measured in seconds rather than months needed for the MD or the detailed Poisson-Nernst-Plank methods (Jakobsson, 1998).

	THE MODEL SYSTEM

TOP ABSTRACT INTRODUCTION THE MODEL SYSTEM RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Experimental input

The experimental signals reconstructed in the present study are the time-resolved fluorescence measurements that were gathered and presented in our previous study of the PhoE channel (Gutman et al., 1992b). The three-dimensional model of PhoE was based on the protein data bank entry 1PHO.

Modeling of the pyranine in the PhoE channel

The PhoE is an anion-selective channel, sensitive to inhibition by polyphosphates. The inhibition is attributed to the clustering of positive charges in the intra-cavity space (Bauer et al., 1988; Samartzidou and Delcour, 1999); thus, we assumed that the pyranine molecule with its three negative charges would interact with the same domain. Upon addition of pyranine to a detergent-stabilized preparation of PhoE, the dye interacts with the protein in a 1:1 stoichiometry, forming a stable complex with G = 9.8 kcal/mol (Gutman et al., 1992b). The complex was noted to be destabilized by replacing the supporting detergent by an anionic one. Furthermore, at 100 mM NaCl the complex dissociated. On the basis of these observations, we had concluded that the pyranine is held in the channel by electrostatic interactions.

Examination of the charge distribution inside the channel indicates clustering of a positive residue near the eyelet of the channel, suggesting that the pyranine (Z = 3) will favor this region. The most probable placement of the pyranine ion in the PhoE intra-cavity was evaluated in two steps. At first, a qualitative search for the positive surface potential was carried out using the GRASP program. Once the binding domain had been identified, a detailed energy minimization was carried out, using the Discover module of the InsightII program.

Fig. 1 depicts the surface potential calculated for the inner surface of the PhoE channel. Most of the surface is colored in red, corresponding with the negative surface charge. The positive regions in the channel (in blue) are grouped in defined clusters, located close to the L3 loop that protrudes into the lumen. The distance between positive domains is compatible with the inter-sulfonate distance of the pyranine molecule.

View larger version (152K): [in this window] [in a new window]   FIGURE 1   A GRASP presentation of the surface electrostatic potential of PhoE (left panel) and of the protein-pyranine complex (right panel). The pyranine, colored in yellow, is located in its most probable location. The upper panels (left and right) depict the electrostatic potential on the protein's surface looking down the channel from the extracellular side. The other panels depict the channel when cut in half through the YZ plane, as marked by the line in the top image. The color scale of the potentials is given at the bottom of the figure.

The pyranine's binding site was calculated by the Discover module of the InsightII program, using the consistent valance force field (CVFF) without any cross-terms. The electrostatic forces were calculated assuming that the dielectric constant of the water inside the cavity was _intra-cavity = 40 (Breed et al., 1996). The pyranine molecule was first inserted in silico inside the PhoE channel, close to the positively charged domain deduced by the GRASP program, and then subjected to structural energy minimization. The first 100 steps were carried out by the steepest descent method, then the PhoE backbone was fixed, and the pyranine molecule was allowed to explore the intra-cavity space using the conjugated mode. The calculations were repeated several times, with varying initial positions and orientations, before the pyranine molecule was allowed to probe the space. In all cases, the final positions of the pyranine molecule were within ± 0.5 Å. Repetition of the calculation with _intra-cavity = 50 (the value that best fitted the experimental observation) did not change the pyranine location in the channel.

Calculation of the electrostatic potential by the DelPhi program

The calculations of the intra-cavity electrostatic potential were carried out by a DelPhi (Honig and Nicholls, 1995) program that was upgraded to account for two dielectric constants assigned for the aqueous phase. The bulk water was set to have  = 80, while the dielectric constant of the aqueous phase in the cavity was an adjustable parameter that varied from _intra-cavity = 2 up to 80. The dielectric constant of the protein was set to be  = 2. The boundary between the protein and the water was defined by rolling a ball r = 1.4 Å over the van der Waals boundary of the protein. The pyranine molecule in its Z = 4 state was placed in the coordinates that were determined as described above and on the assumption of the charge distribution of the ground state molecule.

The calculation of the electrostatic potentials by the DelPhi program utilized the partial charges of the atoms as defined by the PARSE potential set. (Sitkoff et al., 1994, 1996). The size of the PhoE trimer is prohibitively large for mapping its intra-cavity electrostatic potential, which necessitates a very large grid. For this reason, the calculations were carried out for a single monomer, using a field of 169³ grid points. The calculations were initiated at a spacing of 3 Å between two grid points, and the distance between the points was gradually reduced to a full convergence with a resolution of 3 grid points/Å.

The origin of the coordinates (0,0,0) was placed at the oxygen atom of the proton-releasing hydroxyl moiety of the bound pyranine. The mapping was carried out over the whole space of the channel, excluding the space taken by the pyranine molecule plus one solvation layer, i.e., a sphere with a radius of r₀ = 6 Å built around the center of the pyranine molecule.

Karshikoff's calculations (Karshikoff et al., 1994) of the intra-cavity electrostatic potential only considered the fully ionized residues, and the pK values of some groups were modulated to eliminate very steep gradients. The present calculations were carried out at a higher resolution, accounting for the partial charges of all atoms (Sitkoff et al., 1996). Accordingly, all ionizable residues were considered to be in their charged state, in accordance with the standard pK values and the pH of the reaction. It was reasoned that residues of opposite charge, which are at close proximity, will neutralize each other, eliminating the necessity to introduce pK shifts. The ionic strength of the solution was set as 10 µ M, compatible with the experimental conditions used for the laser-induced proton pulse experiments (Shimoni et al., 1993).

The electrostatic potential within the channel space was plotted with a 1 Å resolution and presented as a set of maps representing the electrostatic potential of a 0.3-Å-thick slice perpendicular to the z axis of the channel, which are 1 Å apart. The point with the lowest electrostatic potential was defined for the purpose of the calculation as having E_i = 0 and all other sites were referred to this value.

Defining the trajectory in the channel space

All grid points that were within the boundaries of the channel space were screened, searching at each slice for the grid point with the lowest electrostatic potential (E(i)_min), and their coordinates were recorded.

The most probable trajectory was determined by threading the coordinates of the minima of the consecutive slices, generating an array that extends from the surface of the pyranine molecule up to the opening of the channel to the bulk. On the extracellular section of the channel, the trajectory was initiated at Z = 8 and extended to Z = 22. On the periplasmic side of the channel, the trajectory was initiated at Z = 5 and extended up to Z = 25, where the channel space was regarded as merging with the bulk.

The definition of the transition probability

The probability of a particle stepping between two adjacent loci is a function of the diffusion coefficient of the particle (D_i), the distance between the loci (x), the difference of the electrostatic potential (E) and an entropic term that corresponds to the ratio of the equipotential area of the two loci (S₁/S₂).

(1)

(D_H⁺) in Eq. 1 corresponds with the diffusion coefficient of the proton, which was used as adjustable parameter. x is the width of the increments along the z axis (taken as 1 E) and the entropic terms and gradient of the electrostatic potential are as described above.

Calculation of the intra-cavity transition probabilities

The geminate recombination algorithm of Agmon (Pines et al., 1988) can calculate the propagation of an ion in a defined space under the influence of two forces: the gradient of the electrostatic potential (E/x) and the variation of an entropic term that accounts for the tendency of the particle to sample isopotential sites. The original formalism of Agmon (Agmon et al., 1988; Agmon, 1983; Gutman et al., 1989, 1992a; Pines et al., 1988; Rochel et al., 1990; Szabo et al., 1988), was intended to calculate the electrostatic gradient and entropic terms for a spheric symmetric reaction space. In the present case, the reaction space is not spheric; the electrostatic potentials are affected by the multitude of charges dispersed along the reaction space and the entropic term is a function of both the channel's geometry and the electrostatic potentials, as detailed in Fig. 2.

View larger version (12K): [in this window] [in a new window]   FIGURE 2   A schematic energy profile of successive slices perpendicular to the z axis. (A) The potential of the i slice is lower than that of the next one. As a result, as the proton advances to slice i + 1 it will scan all surface points whose potential is lower than the minimum of the i + 1 slice. The scanned area is marked by hatching. (B) The propagation of the proton along the z axis is downhill. In that case the area scanned by the proton is bound in the range Emin to Emin + kBT.

The figure depicts two scenarios of proton transfer between adjacent slices. In one case (Fig. 2 A), the electrostatic energy of slice (i + 1) is higher than that of the preceding one (E 0), and as the proton propagates toward the bulk, it has to diffuse against the electrostatic potential. As the proton advances by a random Brownian mechanism from slice i to slice i + 1, it will also probe all grid points in slice i having the potential E < E. Thus, the area in slice i that was probed by the proton will vary with E and with the shape and steepness of the electrostatic potential of the i slice. When the energy well is shallow and E is large, the proton's probability density will be smeared over a larger section of slice i, reducing its probability to propagate to the i + 1 position along the z axis. The spreading of the proton's probability density over the i slice is comparable to the entropy term in the transition probability equation. Accordingly, the number of grid points surrounding the minimum (S_i with an area of 1 Å² each) that have an energy in the range E(min)_i  E  E(min)_i+1 were counted and used for the calculation of the transition probability of proton transfer between adjacent slices. Fig. 2 B depicts a situation where the proton propagates down the potential gradient, from slice i toward i + 1, where the energy is lower. In this case, the area probed by the proton was defined as within the range E(min)_i  E  (E(min)_i +1 kT). In all cases where the equipotential area was smaller than the area of a single water molecule (6 Å²), the S_i term was set to be 6 Å².

Propagation of particle along a trajectory

The propagation was carried out by the program of Agmon (Krissinel, 1996) for reconstruction of geminate recombination processes using its FORTRAN (double precision) mode.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THE MODEL SYSTEM RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The most probable location of pyranine

The pyranine molecule, characterized by its molecular coordinates and the partial charges, was placed in silico inside the channel, in the region where the GRASP presentation indicated a positively charged surface. The system was allowed to relax, as described in The Model System, to the most stable configuration. The process was repeated from different initial positions or orientations of the pyranine, and in all cases the position of the dye converged into the same position, as depicted in Fig. 3. The most stable configuration of the pyranine in the PhoE channel is in the eyelet region, with 8 positive (K314, K131, R132, R82, K80, R42, K16 and K18) and 4 negative (E117, D121, D113 and E64) residues within 15 Å from the dye molecule.

View larger version (84K): [in this window] [in a new window]   FIGURE 3   The most probable location of pyranine in the PhoE channel as seen from the periplasmic side of the channel. The protein is schematically represented by the yellow plats for the -sheets, green lines for random coils, and a red cylinder for the  helix. The pyranine is given in atomic detail. The charged residues located within 15 E from the pyranine molecule (blue for basic and red for acidic residues) are presented in ball-and-stick format.

Molecular modeling of the protein-dye complex indicated that the passage through the eyelet was not blocked and the proton released from the excited dye could propagate either to the extracellular or the periplasmic sections of the channel. Thus, the reconstruction of the fluorescence decay curve must account for the possible propagation of the proton along two parallel pathways (extracellular and periplasmic).

The electrostatic potential inside the channel is a function of all charges, including that of the bound pyranine molecule. The pyranine neutralizes some of the positive domains in the eyelet region and, through this, enhances the negative potential of the channel, as evidenced by the reddening of the channel's surface (compare the right and left panels in Fig. 1). Yet, as the GRASP presentation emphasizes the charge density on the surface of the protein and the proton can be dispersed in the total space of the channel, the electrostatic potential in the channel had to be precisely mapped before any reconstruction of the observed signal could be analyzed.

Mapping the intra-cavity electrostatic potential

According to the geminate recombination algorithm (Agmon et al., 1988) the proton is released to the solvent on the surface of a reaction sphere having, in the case of pyranine, a radius of 6 Å. In the present study we adopted this value and considered a proton located within 6 Å from the center of the pyranine molecule to be in the bound state (OH*). The electrostatic potentials were calculated from the surface of the pyranine's reaction sphere to both sections of the ionic channel until it opened to the bulk.

The electrostatic potentials were calculated by the DelPhi program as detailed in The Model System. After the insertion of the pyranine (Z = 3) into the channel, we noticed the presence of three cation traps. These were narrow potential wells, attracting a positive charge by > 10 kT. Two sites were close to the pyranine molecule, while the third one was closer to the vestibule of the periplasmic side. Each of the wells was neutralized by an explicit Na⁺ ion, and the electrostatic potential mapping was repeated.

The electrostatic potentials, as derived by the DelPhi program for Z = 4, are presented as two-dimensional potential maps arranged in sequence along the z axis of the channel (Fig. 4). Each frame represents the value of the electrostatic potential on the xy plane at the z_i coordinate. The calculations were carried out up to a final spatial resolution of 0.33 Å. For presentation, the potential maps are ordered from the innermost one, which is in contact with the first solvation layer of the bound dye, all the way to the bulk. The separation along the z axis is by 2-Å intervals. The protein matrix is marked by dots at 1 Å resolution. The potentials of the grid points in the aqueous phase are color-coded, as shown in the inset to the figure. The left panel depicts the sections along the extracellular section, initiating at Z = 8 up to 22. All along the channel, the maps are characterized by a steep electrostatic gradient that is perpendicular to the z axis. This feature was already noticed by Karshikoff et al. (1994), who only accounted for the net charges of the ionized residues. The steep gradient confines the proton to move along a rather narrow path, adjacent to one side of the channel's wall, indicating that only a small fraction of the intra-cavity space functions as an ion duct (for further discussion, see the accompanying manuscript (Bransburg-Zabary et al., 2002)).

View larger version (179K): [in this window] [in a new window]   FIGURE 4   The maps of the electrostatic potential inside the channel as calculated for the PhoE-pyranine complex. Frames (A) and (B) depict the potentials in the extracellular and periplasmic sections of the channel, respectively. Each image maps the potential in a 0.33 Å thick slice. The slices are 2 Å apart along axis z of the channel. The slice on the upper left corners (Z = 8 and Z = 5) are the closest ones to the pyranine molecule and serve as the sites where the proton is initially released. The protein is represented by the dotted region. The red dots at Z = 5 and Z = 8 represent the Na+ ions that were inserted in the potential traps. The potential scale is given by the color code and the units are kBT.

The high-resolution maps of the electrostatic potential can be used for laborious Brownian motion simulations where a charged particle is allowed to probe the entire space, similar to the calculations of Schirmer and Phale (Phale et al., 2001; Schirmer and Phale, 1999). The alternate procedure adopted in the present study is the first to screen for the grid points that offer the least resistance to proton transfer and consider them to be a one-dimensional propagation trajectory, and then to reconstruct the measured signal as a problem of diffusion in a one-dimensional space. The pathways that the proton will follow through the extracellular and periplasmic sections, derived by screening the potential maps, are represented in Fig. 4.

In the bacterial outer membrane the PhoE protein is assembled in trimers, each having a common extracellular vestibule shared by the three monomers. The common vestibule extends down the channel all the way to Z = 8, yet a proton released in the extracellular section will not spill onto the vestibule space, as the fixed charges on the surface retain the mobile charges close to the surface. The potential map at the Z = 8 reveals a very attractive site, colored in red, that will be spontaneously occupied by the free proton. This lowest potential site was defined as the point of reference and its value was set as E  0. The propagation of the proton from this locus will follow the pathway offering the mildest gradient to progress along the z axis rather than on the xy plane at the Z = 8 slice, where the lateral electrostatic field is as steep as ~20 kT/10 Å. At the xy plane at Z = 12 the electrostatic gradient is much milder, yet the proton will not spill onto the vestibule space (as indicated by the blue line) as the barrier of the electrostatic potential on the path to the vestibule is higher than that of the propagation along the z axis. The leakage to the common vestibule takes place at Z  14.

Panel B of Fig. 4 presents the electrostatic potential along the periplasmic section of the channel. The innermost sections have the shape of a torus, but at Z 15 the edge of the channel opened to the bulk. As in this case of the extracellular channel, the gradient of the electrostatic potential will retain the proton inside the channel with no spilling through the opening to the bulk.

The most probable trajectory

The array of the coordinates of the potential minima in the consecutive slices, each characterized by the electrostatic potential, and the respective S_i values, are sufficient to calculate the respective transition probabilities. This trajectory is a one-dimensional diffusion space, extending from the surface of the pyranine first solvation shell to the edge of the channel.

Fig. 5 depicts the trajectories calculated for the extracellular (left) and the periplasmic sections of the PhoE-pyranine complex. (The calculations were carried out with the three Na⁺ ions filling the sites mentioned above). The pyranine, with a charge of Z = 4, is presented in its most probable position inside the protein (outlined by the -barrel structure). The profile of the electrostatic potentials along the two branches of the trajectory are given at the bottom of the figure, and the lowest point, at Z = 8, was defined as E = 0. The first point on the periplasmic section of the trajectory has a significantly higher potential, ~5 kT above the point of reference.

View larger version (43K): [in this window] [in a new window]   FIGURE 5   The most probable trajectories for proton propagation inside the PhoE channel. The protein is represented by the backbone structure of a single monomer. The pyranine is presented, in a space-filled model, at its most probable location at the center of the protein. The most probable trajectory is composed of two routes, the extracellular (purple) and the periplasmic (orange), by R = 1 Å spheres. The lower section of the figure depicts the electrostatic potentials along the two routes as calculated by the DelPhi program. The lowest potential, at Z = 8, is defined as zero.

The trajectory occupies only a very small fraction of the total channel's volume. It weaves its way through the intra-cavity and in some cases, includes odd, unconnected space elements. The isolated spots correspond with local minima on the pertinent slice, but during the repeated forward and backward stepping, the proton may also probe these "isolated" sites. For this reason, even though they are not on the main path, their entropic contribution to the transition probabilities was not neglected.

The potentials calculated within the channel were affected by the value given to the dielectric constant of the intra-cavity water phase. For low  values the potentials were very negative and the gradients were very steep, although the shape of the trajectory was almost the same. The calculation of the electrostatic potential maps and selection of the trajectories were repeated for varying values of the intra-cavity water phase (2    80), and the transition probabilities along the trajectory were calculated and used for the reconstruction of the measured fluorescence decay curves.

Propagation of proton along the trajectory

The calculation of the propagation of the proton is based on the transition probability and probability density of the proton at each site along the trajectory. The values of the transition probabilities (from slice i to i + 1 and backward) along the trajectory are given in Table 1. Each locus is characterized by its entropic term, the electrostatic potential (calculated for _intra-cavity = 50, a dielectric constant that was found to be the most suitable to reconstruct the observed signal (see below)), and a diffusion coefficient that was set to be identical to that of a proton in bulk water (D_H⁺ = 9.3 105 cm² s¹). As will be discussed below, the diffusion coefficient of the proton could not be determined with high accuracy, yet as the contribution of D_H⁺ to the PT value is a linear function and D_H⁺ is constant along the trajectory, the variations in TP due to this term are not big. The transition probabilities are given in frequency (s1) units and their values imply that a free proton will jump between adjacent loci once within a few picoseconds. The fact that the relaxation of the experimental signal extends over >10 ns is an indication that the proton is propagating, forward and backward, in the reaction space, and many events of geminate recombination, followed by dissociation, take place before the excited dye molecule relaxes to its ground state

Table 1 shows the periplasmic and extracellular sections of the channel. At the end of each section, the proton is irreversibly lost to the bulk. Comparison of the transition probabilities along the two sections of the channel reveals some basic differences between them. The extracellular channel is shorter than the periplasmic one, and through almost half of it (from Z = 8 to Z = 14) the electrostatic potential forms a steep energy barrier that impeded the propagation of the proton to the bulk. In this section of the trajectory, the transition probabilities for advancing toward the bulk are smaller than those for the reverse direction. The periplasmic section of the channel is more favorable for proton propagation and a fast escape to the bulk is expected. Consequently, the reconstruction of the observed fluorescence relaxation is a function of the proton's propagation through the two sections of the channel.

The distribution of the initial proton's probability density between the two sections was assumed to be a function of the electrostatic potential at the innermost loci at the periplasmic and extracellular sections, which are the first to accept the proton as it is ejected from the excited pyranine molecule. The electrostatic potential at Z = 5 is 4.63 kT above that at Z = 8, and according to the Boltzmann distribution function (g_f = exp(E/kT), the discharge of the proton toward the extracellular section will be favored by ~100:1 over the other section of the channel.

Reconstruction of the experimental observation

Fig. 6 depicts the experimental fluorescence relaxation signal as measured for the pyranine-PhoE complex. The fluorescence intensity, at any given time interval, is linearly proportional to the OH* population. Thus, the experimental curves in Fig. 6 are identical to the relaxation of the excited pyranine population. To reconstruct the observed signal, both the initial fluorescence intensity and probability of finding an excited pyranine molecule (P_OH*)_t=0 were normalized to 1. The relaxation of the calculated value (P_OH*)_t and the measured fluorescence intensity were superpositioned. To fit the calculated trace over the experimental one the following parameters were treated as adjustable: 1) the rate of OH* dissociation (_f); 2) the rate constant of proton recombination with the excited pyranine anion (_r); 3) the dielectric constant of the aqueous phase inside the channel space (_intra-cavity); 4) the diffusion coefficient of the proton (D_H⁺); and 5) the rate constant of the reaction of the free proton with the acidic residues on the inner wall of the ionic channel (k_as).

View larger version (23K): [in this window] [in a new window]   FIGURE 6   Reconstruction of the measured fluorescence decay dynamics of the excited pyranine inside the PhoE channel. The main frame depicts the experimental signal and the reconstruction dynamics on a linear scale, to emphasize the quality of the fit in the short time range. The inset presents the data on a logarithmic ordinate, demonstrating that the fit is valid down to 0.1% of the initial amplitude, 20 ns after the initiation of the proton dissociation. The lower frame presents the residue (measured signal minus calculated value) expanded fourfold. The simulation was carried out with the parameters listed in Table 2.

A search over the parameter space yielded the combination of parameters that reproduced the experimental signal with the quality presented in Fig. 6. The main frame depicts the two traces (experimental and calculated) on a linear Y scale, emphasizing the fit during the first nanoseconds after the excitation. The residual values (Y_exp  Y_calc) are given at the top of the main frame and the trend exhibits no systematic deviation. The inset in the figure was drawn on a logarithmic y axis, demonstrating that the predicted function follows the measured signal for almost 20 ns, while the amplitude decreases by more than three orders of magnitude. The values that produce the best fit of the experimental signal are listed in Table 2.

Evaluation of the accuracy of the solution

The large number of adjustable parameters may suggest that more than one combination of parameters can reconstruct the signal. For this reason, it was imperative to evaluate the range of variance for each adjustable parameter. The shape of the best-fitting curve, presented in Fig. 6, was found to be controlled by only four parameters: 1) _f, which is the rate constant of proton ejection from the excited pyranine molecule and is a measure of the chemical activity of the water in the reaction space (Shimoni et al., 1993; Gutman et al., 1982b; Huppert, 1982; Riter et al., 1998); 2) _r, which corresponds with the rate of the back reaction, where a proton in the innermost aqueous layer reprotonates the O* (Cohen, 2001); 3) k_as, which is the rate constant for proton uptake by the carboxylates in the channel; and 4) _cavity, which corresponds with the dielectric constant of the aqueous phase in the channel. The effects of variation of each of these parameters on the reconstructed dynamics are presented in Fig. 7.

View larger version (20K): [in this window] [in a new window]   FIGURE 7   The effect of the adjustable parameters on the reconstructed dynamics. In each frame there are three curves that were calculated with the parameters given in Table 2, except for a single parameter that was set to be either 30% higher or lower than the best-fit value. In (A) the rate constant of proton dissociation (f) was varied. To demonstrate that the deviation from the experimental curve is detected right from the beginning of the reaction, the dynamics are presented in an expanded time scale. Frame (B) demonstrates how the proton recombination rate (r) affects the dynamics. Please notice that most initial dynamics are hardly affected by r. Frame (C) demonstrates how the proton interaction with the carboxylates residues in the channel affects the fluorescence relaxation dynamics.

Panels A-C in Fig. 7 depict the experimental signal with three reconstruction curves; the middle one was calculated with the best-fit value. The other two curves were calculated with the parameter set to be larger or smaller than the best fit by a factor of 30%. It is clear from these figures that deviations of 30% are sufficient to spoil the fit. Above 30%, the distortion of the calculated curve is so severe that modulation of all other parameters cannot restore the fidelity of the reconstruction. Panel A portrays the effect of _f; panel B is for _r, and panel C demonstrates the role of the proton carboxylate reactions taking place inside the channel on the observed dynamics.

The effect of the dielectric constant of the intra-cavity aqueous phase, _cavity, was rather complex. At values smaller than _intra-cavity = 40, the electrostatic forces are hardly screened by the medium and very steep gradients were calculated. Even at the level of double precision routinely used for the calculations, the propagation of the proton in the channel with values of _intra-cavity  40 practically failed, as the numeric solution did not succeed in converging. Calculations that were carried out for the range 45    55 could readily fit the observed signal, and the deviation of the calculated function from the measured one exhibited no systematic deviations. Reconstruction of the experimental curves with the electrostatic potentials calculated for _intra-cavity  55 led to a systematic deviation of the computed function from the measured one, where the predicted curve is consistently above the measured one.

The inability to calculate the dynamics at   40 should lead to a conservative conclusion that the dielectric constant of the intra-cavity water may be smaller than 55, yet the dielectric constant of the intra-cavity water can be corroborated by comparison with other microscopic cavities of comparable size. Recently (Cohen, 2002) we had monitored the dynamics of proton transfer from the excited photoacid molecule (2-hydroxyphenol-6,8-disulfonate) in reversed micelles of varying radii. The analysis was based on the very same formalism used in the present study, where the gradient of the electrostatic potential was calculated by solution of the Poisson-Boltzmann equation for a sphere made of a high dielectric matrix surrounded by a low dielectric continuum (Shimoni et al., 1993). It was noticed that, as the radius of the micelle decreased, so did the dielectric constant of the water inside the micelle. For reversed micelles having a radius of r = 11 Å, comparable to that of the PhoE channel, we obtained a value of  = 60 for the dielectric constant of the water inside the micelle. Combining the two lines of evidence, the value _intra-cavity = 50 ± 5 is an accurate characterization of the partially immobilized water inside the ion-conducting channel of the porin-type proteins.

Interpretation of the kinetic parameters

The parameter controlling the proton dissociation kinetics inside the PhoE channel are given in Table 2, and their values are the basis for the following discussion. Of all terms affecting the reaction space, the properties of the intra-cavity water are the most fascinating ones. The inner space of the PhoE channel is large enough to contain ~650 water molecules, and ~270 are in contact with the van der Waals boundary of the protein that modulates their properties.

Water molecules that are in tight contact with well-solvated ions (Huppert, 1982; Pines, 1989; Riter et al., 1998), or are hydrogen-bonded with the headgroup of phosphoethanolamine (Rand et al., 1988), have a reduced capacity to interact with the charge of the ejected proton. Thus, the rate of proton dissociation from the excited pyranine molecule decreases with the activity of the water in concentrated electrolyte solutions (Gutman et al., 1982a,b; Huppert, 1982). Using the empirical correlation between the rate of the reaction and the activity of water, we estimate the inner space to have a value of a_water = 0.93. It should be stressed that activity is a thermodynamic function. When the reaction space is nonhomogeneous, the measured activity corresponds with the water molecules close to the eyelet, where the dye is located. In the distal domains of the channel, the activity of the water is probably even closer to one. The molecular dynamics calculations of Tieleman and Berendsen (1998) indicated that the water molecules that are in contact with the protein were more polarized than those remote from the interface. This grading of the water molecules could not be refined in the present study, as the observation was too long. As a result, the average properties of the water, which were characterized by a single value for the activity of the water and their dielectric constant, were gauged in the whole space.

The dielectric constant of the intra-cavity space, _intra-cavity = 50 ± 5, is significantly higher than that assigned to the aqueous phase in smaller spaces, increasing from  = 8 and a_water = 0.6, in the heme binding site of apomyoglobin (Shimoni et al., 1993), to  = 30-40 (a_water = ~0.8) in the intralamellar space of multilamellar liposomes. The value we had determined is an experimental confirmation of the theoretical prediction of Sansom (Breed et al., 1996), who deduced that the lower dipolar field of the -barrel (with respect to that of the  helix) would allow a higher level of rotational freedom to the water in the channel. To the best of our knowledge, this is the first experimental determination of a dielectric constant of a microcavity where the detailed structure of the space, at atomic resolution, is included in the considerations.

The special features of a pulse experiment

The measurements carried out in the present study differ from the prevalent steady-state dynamics used for the calculation of a single channel conductance. In the initial state of the present study, the proton binding sites are in equilibrium with each other and their state of dissociation is in accord with their pK values and the pH of the solution. The laser pulse put the system into a temporary state of disequilibrium by releasing a single proton inside the channel's space. The introduction of one excessive proton into a space comparable in size to the PhoE channel will lower the formal pH to ~1, well below the pK of any present carboxylate. The encounter of the excessive proton with any of the carboxylates will lead to the formation of the undissociated state of the residues, which will redissociate within a time frame that is much longer than the present observation time. The dissociation time of a proton from an acidic compound is a function of the pK of the residue (Gutman and Nachliel, 1990); for the excited pyranine molecule the reaction is fast, with  ~ 100 ps. For a carboxylate residue having a pK value of ~4-5, the dwell time of the proton will be a few microseconds. This time frame is orders of magnitude longer than the 20-ns observation time. Accordingly, the binding of the proton to any of the intra-cavity carboxylates will eliminate the probability that the O* molecule will be reprotonated before relaxing to its ground state. Once this reaction is considered, the states of the proton in the present system can be classified into four populations: a proton bound to the excited pyranine molecule (OH*), a free proton with the channel's space, a proton bound to one of the carboxylates, and a proton that diffused out of the channel. The temporal distribution of the proton into the four populations is directly calculated by the propagating program, and the results are given in Table 3. To simplify the calculations, the proton binding sites on the channel's surface were considered to be homogeneously smeared over the surface, with no attempt to identify the proton binding sites with specific residues at given coordinates.

Immediately after the excitation of the pyranine, at t = 0, the proton's probability density is at the OH* molecule. At 1 ns after excitation the pyranine is almost deprotonated (only 6.6% of the proton's population is still in the OH* state), while ~60% of the proton's probability density is in the free form. With time, the fraction of free proton decreases rapidly, with a parallel increase in the bound proton's fraction. Finally, at 15 ns, when the process is over, most of the proton's probability density had accumulated in the carboxylate's bound form and only 1% diffused to the bulk. The efficient removal of the free proton population by the reaction with the carboxylates minimizes the effect of the diffusion coefficient on the reconstruction of the observed dynamics. For this reason, the value of D_H⁺ given in Table 2 is set within a wide range of certainty.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION THE MODEL SYSTEM RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

In the present study we reversed the hierarchical order of the random walk simulations. Instead of allowing the particle to wander at random through the whole diffusion space, followed by a statistical screening for the most probable pathway (Phale et al., 2001; Schirmer and Phale, 1999; Crozier et al., 2001; Jakobsson, 1998), we introduced a selective filter that applied the Boltzmann distribution rule to project the most probable trajectory through which the proton propagates. This algorithm proved itself to be capable of reconstructing the complex experimental signal and in parallel, confirmed the prediction of Sansom (Breed et al., 1996) for the dielectric constant of the intra-cavity water. Such a procedure might be regarded as a refinement that is pertinent only to the specific system under study. To evaluate the general applicability of the procedure one has to test the algorithm under more general conditions, replacing the probe particle (a proton in the present case) by either a positive or negative charge and testing the validity of the prediction on other large-pore channels. As demonstrated in the accompanying publication, the PSST algorithm proved itself capable of predicting the single channel conductance of PhoE and other large-pore channel proteins.