Odorant Binding and Conformational Dynamics in the Odorant-binding Protein*

In mammals, the olfactory epithelium secretes odorant-binding proteins (OBPs), which are lipocalins found freely dissolved in the mucus layer protecting the olfactory neurons. OBPs may act as passive transporters of predominantly hydrophobic odorant molecules across the aqueous mucus layer, or they may play a more active role in which the olfactory neuronal receptor recognizes the OBP-ligand complex. To better understand the molecular events accompanying the initial steps in the olfaction process, we have performed molecular dynamics studies of rat and pig OBPs with the odorant molecule thymol. These calculations provide an atomic level description of conformational changes and pathway intermediates that remain difficult to study directly. A series of eight independent molecular dynamics trajectories of rat OBP permitted the observation of a consensus pathway for ligand unbinding and the calculation of the potential of mean force (PMF) along this path. Titration microcalorimetry confirmed the specific binding of thymol to this protein with a strong hydrophobic component. In both rat and pig OBPs we observed lipocalin strand pair opening in the presence of ligand, consistent with potential roles of these proteins in olfactive receptor recognition.

In mammals, the olfactory epithelium secretes odorantbinding proteins (OBPs), which are lipocalins found freely dissolved in the mucus layer protecting the olfactory neurons. OBPs may act as passive transporters of predominantly hydrophobic odorant molecules across the aqueous mucus layer, or they may play a more active role in which the olfactory neuronal receptor recognizes the OBP-ligand complex. To better understand the molecular events accompanying the initial steps in the olfaction process, we have performed molecular dynamics studies of rat and pig OBPs with the odorant molecule thymol. These calculations provide an atomic level description of conformational changes and pathway intermediates that remain difficult to study directly. A series of eight independent molecular dynamics trajectories of rat OBP permitted the observation of a consensus pathway for ligand unbinding and the calculation of the potential of mean force (PMF) along this path. Titration microcalorimetry confirmed the specific binding of thymol to this protein with a strong hydrophobic component. In both rat and pig OBPs we observed lipocalin strand pair opening in the presence of ligand, consistent with potential roles of these proteins in olfactive receptor recognition.
Odorant-binding proteins (OBPs) 3 have been suggested to act as simple passive transporters of volatile, hydrophobic ligands across the aqueous mucus layer, but may play a more active role in which the neuronal receptor recognizes the OBPligand complex. They may also be involved in the deactivation of odorants (1,2). Only a few OBP subtypes have been found in any given organism, indicating relatively low ligand specificity, with dissociation constants in the micromolar range (3,4). The all ␣-helical OBPs of insects have been suggested to possess intrinsic conformational flexibility (for a review, see Ref. 5) that would allow them to modify their conformation on ligand binding, and which would thereby constitute a mechanism for specific olfactory receptor recognition of the ligand-bound pro-tein. For example, an insect PBP-pheromone complex has been suggested from electrophysiological evidence to be responsible for activation of pheromone receptors (6).
Mammalian OBPs fold into an antiparallel ␤-barrel and are members of the lipocalin superfamily (for a review, see Ref. 7). In many studies, this fold has been assumed to be unlikely a priori to undergo conformational modifications. However, the first structure of a mammalian OBP (8,9), that of the cow, revealed a domain-swapped dimer, in which the helix near the C-terminal region of each monomer packed against the ␤-barrel of the other. Subsequent studies of OBPs from other species and mutant proteins demonstrated that the domain swapping seen in bovine OBP is an exception, arising from a different linker length joining the helix to the barrel domain as well as the absence of the C-terminal disulfide conserved in other lipocalins (10,11). Nevertheless, the occurrence of domain swapping itself is a reminder that the lipocalin structure is not inert.
Indeed, other conformational changes in lipocalins appear to occur under certain conditions, and may be involved in signaling the presence of ligand. In rat OBP, changes in measured amino acid pK a values and circular dichroism spectra suggest that conformational change at some level occurs upon odorant binding (12). Larger modifications of the lipocalin framework are not without precedent as well. The ligand-bound form of the retinol-binding protein is recognized preferentially by transthyretin because of a subtle conformational change induced by ligand binding (13,14). However, recent work on this lipocalin has suggested that it may also undergo conformational changes to release the ligand at the cell surface (15). ␤-Lactoglobulin, despite being the prototypical member of the lipocalin superfamily exhibiting the characteristic ␤-barrel structure, nevertheless has a high ␣-helix propensity; in denaturation/renaturation experiments, structures with a higher helical content than the native structure appear as intermediates (16). NMR data from refolding (17) and pressure-induced unfolding (18) studies of ␤-lactoglobulin have also indicated the presence of intermediate structured forms.
The dynamics of OBP and their potential modifications upon ligand binding have been little studied, but it should be kept in mind that any modulation of pre-existing conformational equilibria of the lipocalin framework by ligand binding would confer upon OBPs a more complex role than that implied in a simple passive transport or scavenging model, and could be linked to interactions of the protein with other partners. Data concerning putative interactions between OBP and olfactory receptors are also scarce, but Matarazzo et al. (19), using purified, radiolabeled protein, have presented evidence of selective, nanomo-lar binding of unligated porcine OBP to a human olfactory receptor, suggesting that OBPs may well be involved directly in the signal transduction mechanism. In such a case, ligandlinked conformational changes would have important consequences in the molecular events accompanying the initial steps of olfaction.
We report here molecular dynamics studies of the dynamics of the OBP-odorant binding interaction for rat and porcine OBP. Such calculations provide atomic-level descriptions of conformational changes and pathway intermediates that are difficult to study directly. The odorant molecule thymol was chosen for this study because of its small size and limited flexibility as well as the availability of the crystal structure of its complex with porcine OBP (3). From a series of independent molecular dynamics trajectories, we observed a consensus pathway for ligand unbinding for rat OBP-1F, and calculated the PMF along this path. Titration microcalorimetry confirmed specific binding of thymol to rat OBP and permitted experimental measurement of the binding affinity, which is interpreted in terms of the PMF and the dynamics results. Marked changes in the dynamic properties of OBP are seen in the presence of ligand, and suggest a promising direction for study in understanding the role of OBPs in olfaction.

MATERIALS AND METHODS
Starting Structures-Two OBP protein structures were used in this work: the crystal structure of the porcine OBP-thymol complex (Ref. 3, PDB identifier 1E06) and the rat OBP-1F homology-modeled structure (12) in which hamster aphrodisin (1EFP; 43% sequence identity) was used as template after secondary structure-weighted multiple alignment of six related lipocalins of known structure. The homology model was energy-minimized and verified using Procheck (20,12). Additional ProsaII analyses (21) revealed the entire structure to be in the favored negative-energy region. Further verifications of the homology-modeled structure are described under "Results." Initial thymol placement in the rat structure was obtained by superimposition of the porcine OBP-thymol complex, followed by energy minimization using CHARMm (22).
Force Field and Parameterization-Molecular dynamics trajectories for the rat and porcine OBP-thymol complexes in explicit solvent were calculated using parameter set 27 of the CHARMm package (22) as described below. InsightII (Accelrys, Inc.) was used to calculate the atomic partial charges on the thymol molecule using MOPACK-AM1. As has been pointed out elsewhere (e.g. Refs. 23 and 15), sampling in explicit solvent simulations is generally much more computationally intensive than in implicit solvent approaches. In particular, equilibration of explicit solvent molecules during ligand passage through a deep binding channel, as is the case in the lipocalin systems, depends on highly stochastic events and poses a particular challenge to canonical or microcanonical explicit solvent approaches (see Refs. 24 and 25). Implicit solvent models not only mitigate such solvent sampling problems but also globally increase conformational sampling. Thus, for unbinding trajectories (described below), we used the CHARMm implementation of the SASA force field (26), providing a surface area-based estimation of hydrophobic surface solvation, while retaining the dielectric model and reduced charge approach of EEF1 (27). Use of this force field provides better qualititative matches of the ␣-carbon fluctuation profiles obtained from crystal structure temperature factors than those obtained with simpler approaches (26). Here, fluctuation profiles were seen to be very similar to those obtained from crystallographic data for aphrodisin and porcine OBP (data not shown).
Design of the Restraint Potential-For unbinding studies an umbrella sampling approach was taken, in which a restraining potential V r was added to the internal energy, defined in terms of the distance d of the thymol center of mass (CM) from that of a fixed set of atoms defining the binding site in Equation 1, in which d 0 is a chosen target distance and the force constant k r ϭ 10 kcal Å 2 .
The set of atoms used to determine the binding site was obtained operationally. The average position of the ligand CM was first calculated in an unrestrained 200-ps preliminary trajectory of the rat OBP-thymol complex. The binding site was then defined as the set of atoms contained within a sphere of optimal radius (12 Å) centered at this position in the protein structure obtained at the end of the preliminary trajectory. The optimal radius was chosen empirically to be large enough to minimize fluctuations of the center of mass of the included atoms while not being so large as to extend beyond the protein surface. The resulting restraint potential is thus defined in terms of an internal coordinate, eliminating the need to apply rotational and translational constraints to the protein-ligand system.
Molecular Dynamics Trajectories-CHARMm was used for all molecular dynamics calculations and subsequent structural analyses. In all cases the system (protein alone or protein plus solvent) was initially energy minimized using harmonic restraints about the starting protein atom positions using a force constant starting at 250 kcal/mol-Å, which was reduced by half in a cyclic fashion after each round of minimization until it fell below 10 kcal/mol-Å and removed completely before final minimization. The system was then heated to 300 K in 25-degree steps using a different random seed for each trajectory and equilibrated for 100 ps.
Unrestrained explicit solvent calculations (1-2 ns each, ϳ23,000 atoms depending on the system) were run in the NPT ensemble (1 atm, 300 K) using periodic boundary conditions and an orthorhombic geometry, with a 9-Å nonbonded cutoff and shift electrostatics in the cutoff region (28,29). SHAKE was used to constrain heavy atom hydrogen covalent bonds. Both 1and 2-fs integration time steps were tested in the explicit solvent calculations with no appreciable differences in the stability of the simulations. Each explicit solvent trajectory was calculated using a distinct distribution of counterions (calculated using Solvate, Ref. 30) and random seed. As to the restrained MD used to explore thymol unbinding, each unbinding run consisted of a single, long trajectory (5-15 ns) using a 1-fs timestep and the method of Berendsen et al. (31) for temperature control. The ligand site distance restraint d 0 was maintained at 0 Å throughout minimization, heating, equilibration, and the first 100 ps of production, and was then increased by 0.2 Å every 100 ps. After each modification of the distance restraint, dynamics runs were restarted with the coordinates and velocities preserved from the previous step. The resulting trajectories were quasi-continuous, exhibiting no abrupt changes in system properties upon application of successive umbrella potentials.
Hydrogen Bonding Calculations-Hydrogen bonding between the ligand and the protein were evaluated using a heavy atom distance cutoff of 3.2 Å, an angular cutoff of 45 degrees from linear for the D-H.A atoms and a lifetime detection threshold of 1 ps.
PMF Calculation-The potential of mean force as a function of the ligand site distance was calculated using the umbrella sampling technique coupled with application of the weighted histogram analysis method (WHAM; Refs. 32 and 33), as implemented in a program by Alan Grossfield. A compact description of PMF methodologies is given by Kosztin et al. (25); their use in particular with implicit solvation treatments are shown in ref 34, in which PMFs obtained using both implicit and explicit solvent representations were compared for the description of ion-pair formation. Woo and Roux (35) and Roux and Simonson (36) describe related approaches including the calculation of protein-ligand affinities from molecular dynamic trajectories using implicit solvent representations.
Protein concentration, as determined by UV spectroscopy using a extinction coefficient of 14,173 M Ϫ1 cm Ϫ1 at 276.3 nm, was 20 M in 50 mM phosphate buffer at pH 7.0 and 30°C in the cell. Odorant solution in the syringe (200 M in MeOH 0.2%) was injected in 40 successive 5-l aliquots at 4-min intervals. Data, corrected for heat of ligand dilution, were analyzed as described previously (12).

RESULTS
The rat OBP-1F system ( Fig. 1) was chosen for study of the unbinding process because of the availability of thermodynamic and biochemical characterizations of this system together with a well documented homology-modeled structure (Refs. 37, 38, and 12 and references therein). We performed additional verifications of the homology model by comparing its dynamic properties to those of the homologous porcine OBP as determined by x-ray crystallography (3). For each of the rat and porcine OBP-thymol complexes, several unrestrained molecular dynamics trajectories in explicit solvent (1-2 ns each) were performed, in which the thymol remained in the calyx throughout. The ␣-carbon root mean square (rms) deviation of the rat OBP ␤-barrel with respect to the starting structure was small throughout the trajectories (0.8 -1.8 Å), and very similar to those calculated for the corresponding trajectories obtained with the porcine OBP structure (0.7-1.5 Å). Further, backbone hydrogen bonding patterns obtained from the explicit solvent trajectories for both rat and porcine structures presented very similar characteristics overall (presented below), as did backbone water bridging interaction patterns (data not shown).
Unbinding Pathways-Several types of molecular dynamics trajectories were performed to study the thymol unbinding process in rat OBP. In preliminary runs, a 200-ps unrestrained simulation was performed and used to define the zero position of the restraint potential as described under "Materials and Methods." Then eight individual trajectories (up to 15 ns), each incorporating a series of incrementally increasing distance restraints, were performed to characterize the ligand unbinding path. Finally, a control run was carried out with the protein alone (10 ns). An implicit solvation approach was used for all unbinding studies in order to increase conformational sampling during the unbinding process. The lipocalin barrel remained stable throughout all the trajectories, with only one unbinding run exhibiting ␣-carbon rms deviations greater than 2.5 Å.
In the absence of any angular restraint, the thymol was free to choose the path taken in exiting the protein in each individual trajectory. In all cases the ligand emerged toward the wider end of the barrel, except trajectory 4, in which it traveled in the opposite direction. Pairwise averages of the differences in the paths followed by the ligand center of mass showed that three trajectories were very similar (average distance, Ͻ2.2 Å), and five were within 3 Å of each other. These five trajectories will be referred to collectively as the consensus pathway, which we describe here in detail.
The average position of the thymol center of mass in the different trajectories is shown in Fig. 1, as well as the consensus path. The variation in color from blue to red indicates the relative degree of fluctuation of the thymol center of mass at a given restraint distance in a given trajectory. Fig. 2 depicts the amino acids interacting with the thymol along this path. In the top panel, the principal hydrogen bonding partners are shown. In the lower panel, the interaction surface is shown, with the color code indicating the degree of hydrophobic contact, measured by the number of molecular dynamics steps for which the contact was made. Residues interacting the most strongly with the thymol in the binding region are shown explicitly in this figure. The side chains of Ala 37 and phenylalanines 55 and 88 made extensive hydrophobic contacts with the ligand in the binding region. Asn 86 , a hydrogen-bonding partner farther along the unbinding pathway, also provided hydrophobic contacts in this region. Fig. 3 shows a histogram of the percentage of molecular dynamics integration steps spent in hydrogen bonds at different distances along the unbinding pathway. Little hydrogen bonding occurred at restraint distances smaller than 2 Å, corresponding to the minimum energy position (see PMF discussion below). Indeed, the thymol is almost completely buried in this region (see curve in Fig. 3), suggesting a principally hydrophobic character for the binding interaction. The most extensive hydrogen bonds were formed at larger distances along the unbinding pathway, primarily with the carbonyl oxygens from Leu 35 mentioned earlier and Tyr 82 , as well as with the side chains of Asn 86 , Ser 83 , and Asn 104 (Fig. 2). Many of the observed hydrogen bonding partners form a neck at the entrance to the internal cavity, the base of which narrows to form the actual binding pocket for the aromatic ligand. The shape of the binding pocket is also reflected in the values of the thymol solventaccessible surface measured along the consensus unbinding pathway. As shown by the curve in Fig. 3, there was virtually no exposed thymol surface in the depths of the binding pocket. About 2 Å into the unbinding pathway, however, the accessible surface began to gradually increase up to restraint distances of about 5 Å. The value then remained constant or even decreased with increasing restraint distances until about 6 Å, before increasing again steadily until the thymol exited the protein.
For the large part, the normal to the plane of the thymol aromatic ring was oriented away from the ␤-barrel axis throughout the unbinding pathway (supplementary Fig. S1). This effect stems from the somewhat flattened shape of the lipocalin ␤-barrel, which tends to align the plane of the ring parallel to the axis. No rotational orientation in the plane of the ring was seen, judging from the orientation of the thymol C-O bond throughout the unbinding runs. The absence of a marked orientational preference is consistent with the elevated B-factors seen for different ligands in crystal structures of complexes obtained with porcine OBP (3).
Free Energy Profiles and Affinity-The conformational sampling obtained from molecular dynamics trajectories using umbrella sampling in either implicit or explicit solvent can be used to calculate a free energy profile for a given molecular process, in this case ligand unbinding. This potential of mean force for the unbinding trajectories was calculated with the WHAM procedure, using the overlapping histograms of ligand-site distances obtained from the umbrella sampling procedure. By definition, the PMF must approach zero at a sufficiently large distance separating the ligand from the protein, but inadequate sampling prevents this limit from being attained with the present approach.
The PMF calculated for the pooled trajectories of the consensus path is shown as the heavy line in Fig. 4, after being set to zero at the plateau that was reached at about 12 Å in all the consensus simulations. At this point the ligand had exited the cavity but remained in contact with the protein. Also appearing in this figure are the PMFs calculated individually for the consensus trajectories. Although the overall profiles are related, the PMF from any single trajectory clearly has much larger uncertainty because of the smaller sampling period. The minimum in the PMF was observed at d ϭ 1.2 Å. This free energy minimum does not coincide with the starting ligand position (d ϭ 0), as a result of the fact that the zero position was chosen in a nonphysical way, being assigned by averaging over ligand positions in a preliminary MD run (see "Materials and Methods"). As seen above, the minimum corresponds to a region of the path in which there are extensive hydrophobic contacts (Fig. 2) and very little solvent exposure of the ligand (Fig. 3). At a distance of about 4 Å the ligand appears to be partly stabilized, despite the gradually decreasing packing interactions at that distance. Increasing hydrogen bonding interactions described above may partially compensate the poorer packing around the thymol in this region. Beyond this distance no net stabilizing interactions between the OBP and the thymol ligand were seen.
The results presented up to this point are independent of the absolute value of the unbound end point of the PMF. Knowledge of this value is of interest, however, because it allows one to calculate the binding affinity theoretically (e.g. Ref. 35). In the current study such a calculation cannot be made rigorously, for we cannot attribute the end point value directly, and, further, no quantitative assessment of free energies obtained with the SASA approach has yet appeared. With this caveat, however, it is of interest to note that in Fig. 4 the PMF appears to reach a plateau at ligand site distances of 12 Å and can be reasoned to be essentially zero. This implies that no net change in free energy accompanies further separation of the ligand from the protein, by the following argument. First, the unfavorable free-energy change associated with solvating the remaining 450 Å 2 of hydrophobic surface (protein plus ligand) still buried at the protein-ligand interface can be calculated from the SASA parameters, yielding a value of 5.4 kcal/mol. This will be offset by the favorable contribution arising from the increase in orientational and positional entropy on ligand release, which should be of approximately the same magnitude; for example, the orientational and translational entropic contribution for the binding of benzene to lysozyme has been calculated at 7 kcal/mol (Ref.  39, see also Ref. 40). Thus, setting the PMF to zero at the plateau amounts to assuming a compensation at longer distances between the loss of the remaining hydrophobic contact and the gain of full rotational and translational freedom of the ligand. With this rationale, the PMF can formally be integrating using the standard relation (e.g. Ref. 41) to obtain a value of 1.6 ϫ 10 9 Å 3 . After being related to standard concentrations (42,35) this value yields a dissociation constant of 0.2 M, corresponding to an association standard free energy change of Ϫ9.2 kcal/mol. Whereas clearly dependent on a certain number of approximations, this value falls nevertheless in the micromolar range seen for a variety of related odorants for OBP-1F (12) and porcine OBP (3). More sophisticated methods such as those employed by Woo and Roux (35), involving explicit solvent molecular dynamics calculations and additional restraint potentials, could provide a more reliable value for comparison to experiment, but would require a computational effort that is far beyond the scope of the present study.
We used isothermal titration calorimetry to confirm the specific binding of thymol to OBP-1F (Fig. 4, bottom). The measured binding affinity was 2.3 Ϯ 0.2 M, giving an association free energy change of Ϫ7.8 kcal/mol, comparable to the estimate given above. The corresponding enthalpy and entropy changes were Ϫ5.4 kcal/mol and 8.0 cal/mol-deg, respectively. Compared with other ligands (12), thymol binding is accompanied by a favorable entropy change, indicating a substantial hydrophobic component to the binding. This is in agreement with the results presented above.
Barrel Splitting and Interstrand Hydrogen Bonding-As mentioned earlier, of the eight unbinding trajectories, number 4 was anomalous in that the ligand traveled toward the small end of the ␤-barrel. A critical point in this pathway was the transient opening of a gap between parts of ␤-barrel strands D and E while the rest of the ␤-barrel hydrogen bonding interactions remained intact. Indeed, in this trajectory the minimum distance between the backbone atoms of strands D and E (the amide nitrogens and the ␣ and carbonyl carbons) widened from 4 Å to a maximum of 13 Å, which allowed passage of the thymol beyond residues forming the base of the binding pocket.
Such barrel splitting might be thought to be a rare event, limited to the rather non-representative pathway sampled in trajectory 4. This was not the case. Indeed, in all other trajectories the unbinding runs were marked by transient opening of similar gaps between ␤-strands D and E of the lipocalin barrel (see inset in Fig. 5). The histogram in Fig. 5 shows the distribution of D-E interstrand distances for the consensus path trajectories. A large percentage of the trajectories can be seen to involve significant strand separation, similar to that seen in trajectory 4. The minimum distance between strands D and E either fluctuated around the starting value (trajectories 2 and 7) or made excursions to values of about 7 Å (trajectories 3, 5, 6, 8) or even 9 Å (trajectory 1) before spontaneously returning to the starting value. The opening was transient in all cases, attesting to the robustness of the lipocalin structure. On the other hand, a 10-ns control trajectory of the OBP protein alone showed virtually no strand opening (see Fig. 5), indicating that barrel opening was associated with the presence of ligand.
To try to better understand this phenomenon we examined the hydrogen bonding between pairs of ␤-strands in the lipocalin barrel. The number of these bonds fluctuates throughout a given dynamics run, but the average number of H-bonds for the different strand pairs is consistent across the trajectories, as shown in supplementary Fig. S2. The smallest numbers of H-bonds were formed between strand E with either strand D or F to either side. The DE pair in particular exhibited marked fluctuation in the number of bonds, with half of the trajectories showing an average number close to zero while in the others this number was closer to 1. This is consistent with analyses of the lipocalin structure (43) in which it was noted that strand E contains only a few residues and the angle between strands D and E is nearly 90 degrees.
Transient destabilization of this side of the lipocalin ␤-barrel is consistent with NMR results obtained by the Goto group for ␤-lactoglobulin, using stopped-flow refolding (17) and equilibrium pressure-induced unfolding (18). In both these studies, what was termed the non-core side of the barrel (strands B-E) was seen to be both kinetically and thermodynamically less stable than the core side (strands F-H). We observed a significant anticorrelation (Pearson's r ϭ Ϫ0.38, p Յ 0.05) between the number of H-bonds formed among pairs of core (F-H) and non-core (B-E) strands (supplementary Fig. S3). Such an anticorrelation suggests dynamic compensation between strand interactions around the lipocalin barrel, which would help offset the energetic cost of transient barrel deformations such as that seen above.
Environments of Tyrosines 20 and 78-Spectrophotometric titration experiments (12) showed that binding of linalool to OBP-1F was accompanied by an increase in the pK a values of two tyrosine residues. Based on the static model structure, it was suggested in that work that the two tyrosines affected by binding should be the conserved Tyr 78 and the less conserved Tyr 20 . The present study supports this assignment, revealing the environments of these two side chains to be modified significantly in the unbinding trajectories versus those of the unligated protein. This can be seen in Fig. 6, which shows histograms of the number of hydrophobic atom centers within 6 Å of the tyrosine hydroxyl group in the different trajectories for these two residues as well as for residue Tyr 82 located near the consensus ligand binding pathway but at the surface of the protein. For Tyr 20 and Tyr 78 , the number of hydrophobic atom neighbors increased significantly with ligand present. The average degree of hydrogen bonding of these two residue side chains with protein neighbors showed a corresponding increase (data not shown), coherent with the less-polar environment. In the case of Tyr 78 , the observed modification may stem from this residue position in strand E in the ␤-barrel, whose packing underwent significant dynamic modifications in all of the unbinding trajectories, as mentioned in the preceding section. The polarity of the environment of Tyr 82 , situated near the entrance to the binding cavity, underwent far smaller modifications. A role for this relatively conserved residue in binding is of course not excluded; it simply appears not to contribute to the pK a shifts seen in rat OBP-1F.
OBP-Thymol Dynamics in Explicit Solvent Calculations-One of the motivations for the use of the implicit solvent approach employed in the unbinding studies reported here is that conformational sampling is far more extensive than can be attained with explicit solvent methodologies (36,44). The development of implicit solvent approaches is widespread, and has allowed investigators to address larger systems or longer timescale phenomena such as ligand binding and conformational change with increasing confidence (45,23). However, the most detailed view currently attainable of certain molecular processes at the atomic level is undeniably provided by explicit solvent simulations, despite the limitations that such approaches themselves remain approximations (typically using non-polarizable atoms, for example) and that long time scale events can only be observed infrequently or not at all.
We performed unrestrained, explicit solvent molecular dynamics calculations on the rat and porcine OBP-thymol complexes. These trajectories confirmed the relative instability of the ␤-barrel strand pairing in the noncore side of the ␤-barrel in the two systems. The backbone hydrogen-bonding matrices presented in supplementary Fig. S4, calculated for comparable 1-ns trajectories of the two structures, are globally similar for the two proteins, again attesting to the validity of the homology model for rat OBP-1F. All trajectories indicated weaknesses in the strand pairing in this region (residues 71-90). In the energy-minimized structures of both porcine and rat OBPs, residues of strand E formed hydrogen bonds with their canonical partners in strands D and F on either side. However, at 300 K interstrand hydrogen bonding tended to be maintained for either the DE or the EF strand pairs, but seldom both. Indeed, water insertion led to significant partial strand pair opening in this region at several points in the trajectories (Fig. 7), whereas no such opening was observed in other strand pairs of the ␤-barrel. In the rat, reversible disruption of backbone hydrogen bonds occurred in the DE and EF strand pairs, while EF strand pair disruption was principally observed in the porcine system. Although a measure of the significance of these differences is not possible due to the limited time scales accessible in these trajectories, it can be noted that the sequence of the DE ␤-turn (residues 74 -77) is not identical in the two species: EGNT in porcine OBP, EDGR in the rat. The glycine at position 75 in porcine OBP allows an additional mainchain hydrogen bond to form between Glu 74 and Thr 77 at the center of the DE ␤-hairpin. This hydrogen bond was persistent throughout the trajectories calculated for the porcine OBP (data not shown), and may strengthen the DE strand pair compared with the rat.

DISCUSSION
In these studies we identified and characterized a consensus pathway for thymol unbinding from rat OBP-1F. This information was obtained using a molecular dynamics unbinding protocol that involves no presuppositions about the unbinding pathway, similar in spirit to an approach used by Chau (46). Hydrophobic interactions were seen to dominate in the depths of the binding pocket, while hydrogen-bonding interactions helped offset reduced packing interactions at larger distances. Isothermal titration microcalorimetry confirmed the specific binding of thymol to this OBP, and the energetics of binding are in the same range as those obtained for this and similar ligands in various mammalian systems (12,3).
Despite the overall stability of the lipocalin structure seen in this study, the molecular dynamics trajectories revealed significant flexibility of the lipocalin molecule in the presence of ligand. This flexibility is notable especially for the appearance of open barrel conformational states arising from strand pair separation, which appeared using both implicit solvation and explicit solvent trajectories for OBP-thymol complexes. Such open states evoke the folding intermediates seen in NMR studies of ␤-lactoglobulin (17,18). Molecular dynamics studies of the human serum retinolbinding protein, another lipocalin, also indicated extensive lability in the region of the D, E, and F strands, which similarly was suggested to play a potential role in dissociation of the retinol-RBPtransthyretin complex and release of retinol at the cell surface (15). And in a recent molecular dynamics study of porcine OBP (47), it was suggested that a shift in the structure of the EF turn carrying Trp 82 may permit ligand passage between the D and E strands. The conformational dynamics observed in the present study involve modification of a small number of hydrogen bonding interactions between ␤-strand E and neighboring strands, whose transient disruption may be compensated by increased hydrogen bonding elsewhere around the ␤-barrel.
Interestingly, the consensus pathway observed in the unbinding trajectories involved the thymol exiting the calyx by passing between the EF turn and the loop 1 residues (Fig. 1). That is, the observed barrel opening was not involved in direct ligand passage through the breach in the barrel wall; rather, it appears that the observed strand pair openings in the rat and porcine systems are simply a consequence of the presence of the ligand in the lipocalin interior. A shift in the equilibrium population of these strand linkages induced by ligand binding could affect putative interactions of the OBP protein with other partners, e.g. the olfactory receptor, in an allosteric linkage mechanism (48,49).
Plasticity in OBPs is presumably necessary for the capacity of these proteins to bind and recognize several different odorant molecules. An analogous conclusion was reached from analyses of crystal structures of ligand-bound artificial lipocalins generated by random mutagenesis of the lipocalin framework by Skerra and co-workers (50). On the whole, such plasticity should confer robustness to the lipocalin framework. However, robustness does not imply rigidity. Thus, while results of FT-IR spectroscopy of porcine OBP (51) did not reveal conclusive evidence of conformational modification upon binding of odorant molecules, those authors noted that even small ligand-induced changes may be sufficient for signaling in other systems such as the retinol binding protein. Whether or not the ligand-linked alterations of OBP dynamics seen in this study are sufficient to allow recognition of the ligation state of OBPs by the olfactory receptor remains to be seen. In any event, the presence of such conformational modifications provide additional, unforeseen elements which may enrich potential chemoreception mechanisms such as those proposed by Matarazzo et al. (19).