INTRODUCTION

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS THEORY: ENTROPY ESTIMATE FOR... RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Molecular recognition is essential for many biochemical processes and may be realized by highly specific binding of ligands to their receptors. In the immune system, antibodies play a crucial role by recognizing and specifically binding their antigens. Calorimetric, spectroscopic, kinetic, and structural data on antibody/antigen systems have been accumulated (Padlan, 1994). Antibodies are widely used in biotechnology, pharmaceutics (Dickman, 1998), and as catalysts (Kristensen et al., 1998). For a better understanding of the antigen-binding mechanism insight into its molecular dynamics at the atomic level is essential.

Recent advances in single molecule atomic force microscopy (AFM) and other force probe techniques allow one to study mechanical properties of individual molecules such as the binding forces of protein-ligand complexes (Lee et al., 1994b; Florin et al., 1994; Moy et al., 1994), the enforced unfolding of proteins (Rief et al., 1997a; Kellermayer et al., 1997; Tskhovrebova et al., 1997), or the stiffness of other polymers, such as DNA (Lee et al., 1994a) or polysaccharides (Rief et al., 1997b). Also, and motivating the present work, single-molecule AFM studies on a number of antibody/antigen complexes have been carried out (Hinterdorfer et al., 1996; Dammer et al., 1996; Allen et al., 1997; Ros et al., 1998), and more are under way.

Computer simulations of such AFM experiments (Grubmüller et al., 1996; Konrad and Bolonick, 1996; Rief et al., 1997b; Izrailev et al., 1997; Marrink et al., 1998; Lu et al., 1998; Heymann and Grubmüller, 1999a-c) by means of "force probe" molecular dynamics (FPMD) simulations can provide models that explain the measured forces in terms of molecular structures, unbinding reaction pathways, and interatomic interactions. They thereby complement single molecular AFM experiments with microscopic interpretations. Such simulations can and have to be validated by computing unbinding forces from the simulations and comparing these with experiments. Such a comparison is complicated by the fact that AFM experiments are typically carried out on a millisecond time scale or slower, whereas MD simulations are currently limited to nanoseconds. Because the unbinding force generally depends on how fast the unbinding is enforced, care has to be taken when relating simulations to experiments (Grubmüller et al., 1996; Izrailev et al., 1997; Marrink et al., 1998).

In a prior publication we reported the computation of unbinding forces of an antibody/hapten complex from nanosecond FPMD simulations (Heymann and Grubmüller, 1999a). In particular, we studied the enforced unbinding of a spin-labeled dinitrophenyl (DNP-SL) hapten from the AN02 antibody. The AN02 antibody is one of 12 monoclonal anti-dinitrophenyl spin-label antibodies that have been characterized in terms of their cDNA sequences and binding constants (Leahy et al., 1988). The AN02 antibody is particularly well-characterized by NMR measurements (Theriault et al., 1991, 1993; Martinez-Yamout and McConnell, 1994) and picosecond MD simulations (Theriault et al., 1993).

To rescale the unbinding forces computed from our AN02 FPMD simulations from the nanosecond simulation time scale to the experimental millisecond time scale (and longer), a scaling model has been developed and applied (Heymann and Grubmüller, 1999a). This model considers both friction and activated processes and makes use of the measured spontaneous dissociation rate for the AN02/DNP-hapten complex and the unbinding length that was estimated from the simulations. From this approach, we predicted an unbinding force of 60 ± 30 pN for the millisecond time scale. We expect experimental AFM data on this system to be available soon.

In this paper we focus on the microscopic interpretation of the unbinding forces of the AN02/DNP-hapten complex in terms of physical interactions and unbinding pathways. To that end, we consider both the interactions between the ligand and individual AN02 binding pocket residues and more structural aspects of the unbinding pathways. In particular, we focus on a possible structural heterogeneity of the unbinding pathway, which is assumed to cause significant scatter of unbinding forces. Such scatter has already been observed in a number of systems (Florin et al., 1994; Moy et al., 1994; Grubmüller et al., 1996), but was generally attributed to experimental or statistical error rather than structural heterogeneity. For the bound AN02 complex, structural heterogeneity has been observed by NMR (Theriault et al., 1991; Martinez-Yamout and McConnell, 1994). Particularly, the light chain Tyr-31 was suggested to be rather flexible in these studies, and to contribute to induced fit motions upon binding (McConnell and Martinez-Yamout, 1996).

To investigate the microscopic nature of such structural heterogeneity and to assess its relevance for a proper description of the kinetics of binding and unbinding, we shall attempt to quantify the entropic contribution to the free energy barrier of enforced unbinding that results from that heterogeneity. To this end, principal component analyses of the unbinding trajectories will be used.

The importance of individual residues with respect to ligand binding is traditionally judged from site-directed mutagenesis of putative residues by measuring the change in free energy with respect to the wild type. Typically, that change can also be estimated from the x-ray structure by counting hydrogen bonds, salt bridges, or van der Waals contacts between ligand and binding pocket. Aiming at extending this procedure toward forces, we identified critical hydrogen bonds that are expected to significantly reduce the AN02/DNP-hapten binding free energy.

To answer the question of how those crucial residues also affect unbinding forces, we modeled two mutants of AN02 and carried out another series of FPMD simulations for each of the two mutants. In the first the mutant, the light chain Tyr-33 that interacts with the hapten molecule via a strong hydrogen bond was replaced by phenylalanine, which can not form such a hydrogen bond. The binding constant of that mutant is smaller than that of the wild type by a factor of 10 (Martinez-Yamout and McConnell, 1994). Accordingly, we also expected to find lower unbinding forces. In the other mutant, which has not been studied previously, the interaction between hapten and the binding pocket was increased by replacing Ile-96 by lysine. The lysine is expected to form an additional strong hydrogen bond to the hapten molecule, which should give rise to an increased unbinding force for this mutant.

	COMPUTATIONAL METHODS

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS THEORY: ENTROPY ESTIMATE FOR... RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Simulation systems and FPMD simulation methods

As a start, the x-ray structure of the AN02 F_ab fragment complexed with the DNP-SL ligand (Brünger et al., 1991) was taken from the Brookhaven Protein Data Bank (Bernstein et al., 1977), entry 1baf. All molecular dynamics simulations were performed with the parallel MD program EGO (Eichinger et al., 2000) (Eichinger, M., H. Grubmüller, and H. Heller. 1995. User Manual for EGO_VIII, Release 2.0. Electronic access: www.mpibpc.gwdg.de/abteilungen/071/ego.html) and the CHARMM force field (Brooks et al., 1983). EGO allows efficient computation of Coulomb forces using the "fast multiple time step structure adapted multipole method" (Eichinger et al., 1997; Grubmüller and Tavan, 1998), so that no artificial truncation of these forces had to be applied (Brooks et al., 1983). For the DNP-hapten, partial charges were calculated with UNICHEM (Cray Research Inc. 1997. UniChem 3.0. 655 Lone Oak Drive, Eagan, MN 55151.) (see Fig. 1) using MNDO94, the AM1 Hamiltonian, ESP-fitting, and a dielectric constant of one. All force constants were taken from the CHARMM force field (Brooks et al., 1983). All simulations were carried out with an integration step size of 10¹⁵ s. The length of chemical bonds involving hydrogen atoms was fixed (McCammon et al., 1977), and nonpolar hydrogen atoms were treated through compound atoms (Brooks et al., 1983). No explicit term for the hydrogen bond energy was included within the force field.

View larger version (24K): [in this window] [in a new window]   FIGURE 1   Structure of the spin-labeled dinitrophenyl hapten; shown are names and partial charges as computed with UNICHEM and used in the simulations.

A model for the solvated protein-ligand complex was built by surrounding the x-ray structure with a droplet containing 13,461 TIP3 (Jorgensen et al., 1983) water molecules, 34 Na⁺ ions, and 41 Cl ions at physiological concentration using the program SOLVATE (Grubmüller, H. 1996. Solvate: a program to create atomic solvent models. Electronic access: www.mpibpc.gwdg.de/abteilungen/071/hgrub/solvate/docu.html) [see Fig. 2 A]. The ions were placed according to a Debye-Hückel distribution governed by the protein and hapten charges. The excess of seven Cl ions served to compensate the net positive charge of the complex. To keep the number of solvent molecules in the system as small as possible, a nonspherical solvent volume was chosen. The solvent surface was defined such that the minimum distance between the protein surface and the solvent surface was 20 Å near the binding pocket region and 12 Å elsewhere. The obtained simulation system comprised a total of 44,571 atoms.

View larger version (75K): [in this window] [in a new window]   FIGURE 2   (A) Model of the solvated AN02 (red)/DNP hapten (green) complex used in the molecular dynamics simulations described in the text. The complex has been surrounded by water molecules (blue); ions (yellow) were placed according to a Debye-Hückel distribution governed by the protein and hapten charges. (B) General setup of the force probe molecular dynamics (FPMD) simulations. Shown is the AN02 Fab fragment (ribbon plot) and the hapten molecule (ball-and-stick model). The oxygen atom O2 was subjected to a harmonic potential (spring), which was moved in the pulling direction (arrow) with constant pulling velocity vcant.

To counterbalance surface tension and to prevent evaporation of water molecules, all surface water molecules were subjected to a "deformable boundary" (Brooks and Karplus, 1983), which has been adapted to water solvent and approximated by a quartic polynomial (Kossmann, 1997). Additionally, all water molecules at the droplet surface were coupled to a heat bath of 300 K through stochastic forces obeying the fluctuation-dissipation theorem (Reif, 1965). A coupling constant of  = 10 ps¹ was used. All other atoms were weakly coupled to a heat bath ( = 1 ps¹) via velocity rescaling (Berendsen et al., 1984).

After energy-minimizing the whole system by 1000 steepest descent steps, the system was equilibrated for 1500 ps at 300 K. During equilibration, the stability of the complex was monitored via its root mean square (rms) deviation from the x-ray structure ("forward deviation") and, as suggested by Stella and Melchionna (1998), also from the structure after equilibration ("backward deviation"). Following the suggestion of one of the referees, structural motions observed during equilibration were compared to structural flexibility obtained from a comparison of six related x-ray structures with each other using the program DynDom (Hayward and Berendsen, 1998).

For the rms computations four sets of heavy atoms were considered, namely those 1) of the complete protein-ligand complex, 2) of the variable, and 3) of the constant domain regions, and 4) of the binding pocket, respectively. Here, the binding pocket (referred to as "B" below) was defined as those residues that interact with the hapten molecule during the unbinding process, together with those parts of the hypervariable loops that are located close to the hapten molecule. These comprised residues 30-35(L), 47-52(L), and 86-99(L), as well as residues 30-36(H), 48-55(H), and 97-105(H), respectively. Here, and also in the following, (L) and (H) refer to the light and heavy chain.

The subsequent FPMD simulations were carried out as sketched in Fig. 2 B: in close resemblance to the pulling forces exerted in the single-molecule AFM experiment by the cantilever, the spin-labeled oxygen atom O2 of the hapten was subjected to a "spring potential" V_spring,

(1)

acting in a pulling direction on the z-coordinate z_O2 of atom O2. Here, z_cant(t) is the cantilever position. A spring constant of k_cant = 280 pN/Å was used throughout all simulations.

In the figure, V_spring is symbolized by a spring. Note that V_spring as defined above does not restrain sideward motions of the O2 atom, i.e., perpendicular to the z-direction, because such motions are also essentially unconstrained in the experiments. At the beginning of each FPMD simulation, the minimum position, z_cant(0), of V_spring was placed at the position of atom O2, z_O2(0). In the course of each simulation, V_spring was moved in the pulling direction with a constant pulling velocity v_cant (arrow), i.e., z_cant(t) = z_cant(0) + v_cantt.

To avoid larger sideward drift of the hapten, atom O2 was additionally subjected to a weak harmonic potential (spring constant: 14 pN/Å) perpendicular to the pulling direction. Finally, the center of mass of the protein (with hapten excluded) was kept in place by a relatively stiff harmonic potential (force constant: 2800 pN/Å). This setup allowed the protein to adapt to the enforced hapten unbinding, e.g., by rotations or intramolecular conformational motions like induced fit motions, also in close resemblance to typical AFM experiments.

During enforced unbinding, the pulling force F_pull(t) = k_cant[z_cant(t)  z_O2(t)] acting on the oxygen atom was calculated and recorded every 100 fs. Thus, a "force profile" for each FPMD simulation was obtained as a function of cantilever position z_cant(t). An unbinding force was defined as the maximum of the respective force profile, and the position of the O2 atom at that maximum was recorded to obtain an "unbinding length" (Heymann and Grubmüller, 1999a). The relatively stiff spring constant of k_cant = 280 pN/Å for the spring potential was chosen to provide high spatial resolution of the force profile of k_BT/k_cant  0.5 Å while allowing thermal fluctuations of about the same size (Heymann and Grubmüller, 2000). High-frequency fluctuations of atom O2, which were artificially caused by the harmonic spring potential, were filtered by properly smoothing the force profiles (Grubmüller et al., 1996).

To characterize the velocity dependence of the unbinding forces and to allow proper statistical analysis of the unbinding pathways, a total of 58 FPMD simulations were carried out. For computation of the unbinding force as a function of the pulling velocity, 29 simulations were used. For these simulations, pulling velocities in the range between 0.1 and 50 m/s were chosen requiring simulation lengths between 30 and 7000 ps to complete the unbinding process in each case (Heymann and Grubmüller, 1999a). To study the structural heterogeneity of the enforced unbinding pathways, 20 FPMD simulations of 300 ps length each were carried out with identical pulling velocity, v_cant = 5 m/s. These simulations differed in their initial conformations, which were picked randomly from the last 820 ps of the 1500-ps equilibration trajectory.

Based on the equilibrated wild-type structure (taken from the equilibration phase at 1300 ps), two mutants of the AN02/DNP-hapten complex were modeled. Here our goal was to identify two mutation sites with opposite effects: the first one should decrease, the second one should increase the unbinding force. To that end, we first selected a hydrogen bond between Tyr-33(L) and the DNP molecule which, as judged from our prior FPMD simulations, strongly contributed to the unbinding force. That bond was then removed by a Y33F(L) mutation. The second mutant, I96K(L), was obtained by identifying the oxygen atom O25 of the DNP ring as a putative hydrogen bond acceptor, and by changing a suitably located isoleucine residue to lysine as a hydrogen bond donor.

To generate a model for the Y33F(L) mutant, the Tyr-33(L) side-chain hydroxyl group was removed and the force field was changed accordingly using the XPLOR (Brünger, 1988) phenylalanine force field parameters. The system was subsequently equilibrated for 50 ps.

The I96K(L) mutation required particular attention because here a charged NH group had to be inserted. This was achieved using the molecular editor of the software package Quanta (Molecular Simulations Inc., 1997). Care was taken that no overlap with neighbored amino acids occurred. To compensate the positive net charge of lysine, a nearby sodium ion was removed. The resulting structure was equilibrated for 100 ps.

Starting from the respective equilibrated structures, each of the mutants was subjected to 22 FPMD simulations to allow accurate comparison of unbinding forces between the mutants and the wild type. Pulling velocities from 50 m/s to 2 m/s were used.

Analysis of FPMD unbinding simulations

For each of the FPMD simulations, changes of the interaction network between the hapten molecule and the AN02 F_ab fragment during unbinding were recorded and characterized by computing selected observables as a function of the cantilever position z_cant. These were the total electrostatic (Coulomb) and van der Waals (Lennard Jones) interactions 1) between hapten and AN02 and 2) between hapten and selected individual residues of the AN02 binding pocket, and formation and rupture of 3) hydrogen bonds and 4) water bridges between hapten and the relevant AN02 binding pocket residues. As relevant binding pocket residues we selected those residues which, in the course of the FPMD simulations, show significant electrostatic or van der Waals interactions with the hapten. Those were Trp-90(L), Tyr-33(L), Tyr-31(L), Asp-49(L), Ile-96(L), Gln-88(L), Trp-100(H), Tyr-51(H), Pro-101(H), and Tyr-54(H), and will be referred to as "R" below.

From the FPMD trajectories the electrostatic and van der Waals energies, respectively, among all atoms of the DNP molecule and all atoms of the protein and among all atoms of the DNP molecule and all atoms of the selected binding pocket residues were computed using XPLOR (Brünger, 1988). A cutoff of 10 Å was used for these energy computations. Energies were computed once per picosecond by averaging the values of 10 snapshots taken every 100 fs for each picosecond. To assess the integrity of the binding pocket during equilibration, this analysis has also been carried out on the 1500-ps equilibration trajectory for comparison.

Strength, formation, and rupture of hydrogen bonds between hapten and the binding pocket residues in the course of the unbinding process was monitored using a combined geometric and energetic criterion for the identification and quantification of hydrogen bonds (Kitao et al., 1993). Accordingly, a hydrogen bond was identified if the involved heavy atoms were closer to each other than 4 Å and if the sum of van der Waals and electrostatic energy was negative and its absolute value >1 kcal/mol. This sum was used to measure the strength of the identified hydrogen bonds. Water bridges between the hapten molecule and binding pocket residues were identified by requiring a water molecule to simultaneously form two hydrogen bonds, one to the hapten molecule and one to a binding pocket residue.

To structurally characterize the unbinding pathway of the hapten molecule and induced-fit motions of the whole binding pocket "B," FPMD trajectories were analyzed in detail. Particular attention was paid to ligand motions relative to the light chain residues and to the heavy chain residues of the binding pocket (cf. Fig. 13). Only trajectories from FPMD simulations with pulling velocities (v_cant  2 m/s) were considered for this purpose because those frictional effects have been shown to be negligible. These simulations are thus assumed to best describe enforced unbinding as seen in the AFM experiments. Figures were prepared with Molscript (Kraulis, 1991) and MOLMOL (Koradi et al., 1996).

Principal component analysis (PCA) (Mardia et al., 1979) was used to estimate the size of the region in configurational space (of both ligand and binding pocket residues "B") covered by the ensemble of unbinding pathways obtained in our FPMD simulations. Generalizing the entropy estimate of protein structure ensembles suggested by Karplus and Kushick (1981) and improved by Schlitter (1993) to ensembles of trajectories, the PCA results were used to estimate the entropic contribution to the free energy along the unbinding pathway and, in particular, to the free energy barrier that governs the AN02/DNP-hapten binding kinetics. This procedure is described in the subsequent Theory section.

A total of 20 FPMD simulations has been used for the entropy estimate, each of 300 ps length. From each trajectory, 3000 coordinate sets (one per 100 fs) were taken for the PCA. These sets included the whole binding pocket "B" and the hapten molecule, comprising a total of 475 atoms.

As illustrated in Fig. 3 and explained in the Theory section, the coordinate sets were split into n = 131 overlapping subsets according to the reaction coordinate z_cant(t) using intervals I_i(i = 1, ... , n) of 2 Å width each. The intervals were shifted along z_cant(t) in 0.1 Å steps in the range between 0 Å (bound state) and 15 Å (unbound state). The obtained 131 coordinate subsets of 400 coordinate sets each were then evaluated using Eq. 11.

View larger version (29K): [in this window] [in a new window]   FIGURE 3   Sketch of the unbinding pathways of the hapten molecule out of the AN02 binding pocket. The volume covered by all unbinding pathways in a given interval Ii of cantilever positions zcant is used to estimate the entropy of the complex.

The relative contribution of the individual degrees of freedom to the configurational entropy change was calculated from the partial entropy differences S_m(z_cant) with m = 50, 55, 60, ... , 100. Larger values for m were not considered because at the chosen spatial resolution of v_cantt = 2 Å only 400 coordinate sets were available for each z_cant value.

	THEORY: ENTROPY ESTIMATE FOR REACTION PATHWAYS

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS THEORY: ENTROPY ESTIMATE FOR... RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

We assume a given ensemble {x⁽ⁱ⁾(t)}_i=1...n of n MD trajectories along a reaction pathway with a Boltzmann ensemble {x⁽ⁱ⁾(0)}_i=1...n of initial conditions (see Fig. 3). Here, x⁽ⁱ⁾(t)  ^3N is the configuration vector of trajectory i at time t, where N is the number of atoms considered. Using this ensemble of trajectories, and generalizing the usual expression for the entropy in the canonical ensemble (Reif, 1965), we want to estimate how the entropy of this system changes in time as the system is forced to move along a reaction pathway,

(2)

where (x, t) is the time-dependent configuration space density that generates the trajectory ensemble {x⁽ⁱ⁾(t)}. In the above equation, the integral is taken over the full configurational space. (Because velocity-dependent forces are absent, the contribution of the momenta to the free energy is a constant and is thus irrelevant for the present purpose.)

Note that for the case at hand the ensemble moves along the unbinding reaction pathway with constant velocity, z_cant = tv_cant, hence S(t) will be interpreted as the (configurational) entropy contribution to the free energy landscape along the unbinding reaction pathway as defined by the reaction coordinate z_cant. Note also that, in defining S(t) via Eq. 2, we assume that all degrees of freedom perpendicular to the reaction pathway are in equilibrium at all times; the same assumption is made in Kramers' transition state theory (Eyring, 1935; Kramers, 1940).

For the stationary case, entropy estimates are available. As implicitly assumed by Karplus and Kushick (1981), (x) can be approximated by a multivariate Gaussian function (x) centered at x₀ (Grubmüller, 1995),

(3)

with partition function

(4)

Here, A  ^3N×3N is a symmetric, positive semidefinite matrix which defines the shape of the Gaussian. It is derived from the covariance matrix of the coordinate sets (Gardiner, 1985), C,

(5)

and can therefore be computed from an MD trajectory.

Within this approximation the (classical) entropy S = k_B  d^3Nx(x) ln (x) can be computed (up to an additive constant) from the determinant of the covariance matrix (Karplus and Kushick, 1981), S = 1/2k_B ln |C|, or, equivalently, from the product of the (nonzero) eigenvalues of the covariance matrix.

That estimate, however, suffers from instabilities caused by small eigenvalues that are unphysical in that these instabilities would vanish in a quantum-mechanical treatment of the above approximation. A correction has therefore been proposed (Schlitter, 1993),

(6)

which provides a better entropy estimate, and also cures the instabilities. Here, M is the diagonal mass matrix that contains the atomic masses,  is Planck's constant, and e = 2.718 ... is the Euler number.

We shall use Eq. 6 to estimate the time-dependent entropy S(t) defined in Eq. 2. However, a straightforward generalization of Eq. 6 is hampered by the fact that, because of the small number of trajectories usually available (n = 20 in our case) for each instance in time t_j, the covariance matrix computed from {x⁽ⁱ⁾(t_j)} would be inaccurate due to statistical error. Furthermore, the covariance matrix would have too few (n  1) nonvanishing eigenvalues, which means that only n  1 out of the total number of 3N  6 internal degrees of freedom would be considered for the entropy estimatemost likely too few.

We therefore rewrite Eq. 2 as

(7)

where (t) is the Dirac delta function. Assuming that the inner integral in Eq. 7 is a smooth function on a short time scale t, (t) can be replaced by a (finite) square function, and an approximation for the entropy is obtained

(8)

Here (t) is the step function, which is 0 for t  0 and 1 otherwise. With the above condition of smoothness, S(t) can further be approximated by the time-averaged value of the inner integral,

(9)

for which Eq. 6 is applicable with the understanding that now the subensemble of structures {x⁽ⁱ⁾(t_j)}_ttj<t+t is to be used for the computation of the covariance matrix C, which thus becomes a function of t. For sufficiently large t, that ensemble is large enough such that a sufficient number of degrees of freedom enters into the entropy estimate.

Using the subensemble of the first interval as the (arbitrarily chosen) zero point, the entropy change along the unbinding reaction pathway is then readily estimated,

(10)

The resulting entropy difference S(z_cant) serves to quantify the entropic contribution to the free energy caused by the structural heterogeneity of the unbinding pathways. Furthermore, by replacing the two determinants in Eq. 10 by products of the respective eigenvalues _l(t) (which we assume for each t to be ordered by decreasing size),

(11)

the relative contribution of each (collective) degree of freedom l to the entropy change S(z_cant) can readily be seen. Because, typically, this relative contribution is large for those degrees of freedom with large eigenvalues (García, 1992; Kitao et al., 1991; Amadei et al., 1993; Schulze et al., 2000) and vanishes for those with small eigenvalues, partial entropies S_m(z_cant) of the first m terms in Eq. 11 can be used as a check for proper convergence.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS THEORY: ENTROPY ESTIMATE FOR... RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Equilibration

Fig. 4 shows the rms deviation for selected regions of the AN02 antibody during equilibration. From the traces for the constant region (top left), the variable region (bottom left), and the binding pocket "B" (top right), a significant increase from the x-ray structure (t = 0) during the first 100 ps of the equilibration is seen. After this short relaxation period, the region of interest, namely the binding pocket region, shows no further drift from the x-ray structure and stabilizes at the relatively small rms value of 1.4 Å. Small drifts are seen, however, for the constant and variable regions, respectively.

View larger version (24K): [in this window] [in a new window]   FIGURE 4   Root-mean-square (rms) deviation of selected regions of the AN02 antibody with respect to the x-ray structure during equilibration. Considered are the constant region (top left), the variable region (bottom left), and the binding pocket region (top right), respectively. Root-mean-square deviations per residue are shown in the bottom right panel where rms deviations of >3 Å are highlighted in black.

Although the stability of the binding pocket region seems to be established, the slow rms drifts of the constant and the variable regions deserved further study. To that end, Fig. 4 (bottom right) plots the rms deviation per residue; rms deviations of 3 Å are highlighted in black. These regions of high mobility are also highlighted in black in Fig. 5, where the AN02 x-ray structure (left) and the structure after equilibration (right) are shown as ribbon plots. As can be seen, all high-mobility regions are loops at the protein surface. In contrast, the  scaffold, and in particular the binding pocket, appear nearly unchanged.

View larger version (53K): [in this window] [in a new window]   FIGURE 5   Comparison of the AN02/DNP-hapten x-ray structure (left) with a representative snapshot after equilibration (right). Regions with an rms deviation >3 Å are shown in black (cf. Fig. 4, bottom right).

An even larger rms drift shows up in Fig. 6, where the rms deviation of the complete AN02/hapten complex during equilibration is shown (heavy line). Here an rms value of 3.1 Å is reached at the end of the equilibration phase, and no clear evidence for convergence is seen. That deviation is significantly larger than that of the individual chains, as seen in Fig. 4 (left), which suggests that this increased deviation is mainly due to very slow motions of the constant region relative to the variable region.

View larger version (16K): [in this window] [in a new window]   FIGURE 6   Root-mean-square (rms) deviation of the AN02/DNP-hapten complex from the x-ray structure ("forward deviation," heavy line) and from the final, equilibrated structure ("backward deviation," thin line).

In this case, the relatively large drift of 3.1 Å seen for the whole F_ab fragment would likely result from (equilibrium) sampling motions of the domains rather than from (nonequilibrium) relaxation motions. To distinguish between these two cases, it has been suggested (Stella and Melchionna, 1998) to additionally consider the backward deviation (thin line in Fig. 6), which measures the deviation from the equilibrated (final) structure.

In contrast to the forward rms deviation, the backward deviation stays small (below 1.7 Å) between 550 and 1300 ps, and exhibits only slight drift during that period. Note also that the abrupt drop to zero between 1300 and 1500 ps must be attributed to equilibrium fluctuations, because otherwise a corresponding rms change would show up during that period in the forward deviation as well, which is not the case. Therefore, following the argumentation of Stella and Melchionna (1998), we assume that the main relaxation motions take place during the first 550 ps of the equilibration phase and that the period after 550 ps, in particular the obvious rms drift, is characterized by slow equilibrium sampling motions.

Closer inspection of snapshots during equilibration (not shown) supports this view and reveals a hinge bending motion of the variable region domain of the F_ab fragment with respect to the constant region domain around the two loops connecting these two regions. The variable region, which contains the binding pocket, stayed unchanged during equilibration.

This observation is in full agreement with the conformational flexibility of the F_ab fragment as obtained from a selection of AN02 homologs (Brookhaven Protein Data Bank entries 1hin:L, 1ifh:L, 1him:J, 1him:H, 1hil:C, and 1hil:A), which have been identified as "structural neighbors" with SCOP (Murzin et al., 1995). From these structures, and in a similar manner as suggested previously (Van Aalten et al., 1997), a hinge-bending flexibility between the constant and variable domain region was identified with angles between 7 and 16° and a hinge axis similar to the one observed in the MD simulation (data not shown). Averaged over all 15 pairs of the above six structures, the heavy atom mean rms deviation is 1.33 Å for the whole F_ab fragment and 0.98 Å and 1.08 Å for the variable and constant domains, respectively, which also documents the significant contribution of the interdomain flexibility to the overall rms deviations.

The interaction network that stabilizes the hapten molecule within the binding pocket has been characterized from the x-ray structure and other studies (Brünger et al., 1991; Anglister et al., 1987; Leahy et al., 1988; Theriault et al., 1991, 1993; Martinez-Yamout and McConnell, 1994). These interactions, shown in Fig. 7, are unchanged after the equilibration phase. The DNP ring of the hapten molecule (gray) is buried in a tryptophane sandwich formed by Trp-90(L) and Trp-100(H), and stabilized by contacts between the aromatic rings. Several hydrogen bonds further stabilize the complex; the strongest are shown in the figure (dashed lines). Here, the hydrogen bond between Tyr-33(L) and an amino group (N8,H29; cf. Fig. 1) at the linkage of the two hapten rings is the most stable. Also observed are hydrogen bonds to Trp-90(L), Gln-88(L), and Ile-96(L), as well as (not shown) a weak one to Tyr-31(L) and a fluctuating one to Asp-49(L).

View larger version (27K): [in this window] [in a new window]   FIGURE 7   Interaction network between bound DNP-hapten (gray) and AN02 binding pocket after equilibration. Strong hydrogen bonds between the hapten and binding pocket residues are indicated by dotted lines. Residues are numbered as in the x-ray structure (Brünger et al., 1991) and differ by one from the numbering used by Martinez-Yamout and McConnell (1994).

Fig. 8 quantifies the above statement. Shown is the strength of the relevant hydrogen bonds and that of the van der Waals contacts between hapten and the relevant binding pocket residues "R" (encoded via the thickness of the lines) during the equilibration period. Although the strength of the interactions fluctuates, as can be expected from the dynamics of the system, no overall and significant change of the interaction "fingerprint" is observed. Consider, e.g., the hydrogen bond to Tyr-33(L), which weakens at 970 ps, but is completely restored afterward at 1150 ps. With a possible exception of the bond to Asp-49(L), which appears somewhat weaker than in the x-ray structure, the interaction network and, thus, the binding pocket region remain intact after equilibration.

View larger version (30K): [in this window] [in a new window]   FIGURE 8   Sum of electrostatic and van der Waals interactions between DNP-hapten and individual residues of the light chain (top half of the plot) and the heavy chain (bottom half) during equilibration. The thickness of the lines reflects the strength of the interactions.

From the above results we conclude that the AN02/hapten bond is described by our simulations to sufficient accuracy. Thus the obtained equilibrated structure provides an appropriate starting point for the subsequent FPMD simulations, the results of which will be described. We also consider the period between 550 and 1500 ps suitable to provide an ensemble from which independent starting structures for otherwise identical FPMD simulations can be extracted as a check for reproducibility and to improve the statistics of the simulation results.

Enforced unbinding of the AN02 antibody/hapten complex

Two sets of FPMD simulations were carried out. To study the variation of unbinding force on the pulling velocity v_cant and to check to what extent unbinding pathways and force profiles depend on v_cant, 29 simulations with v_cant between 0.1 and 50 m/s were carried out. For structural analysis and computation of force profiles, only the simulations with v_cant < 2 m/s were considered, for which friction was found to be negligible (Heymann and Grubmüller, 1999a). Within that range, as already described (Heymann and Grubmüller, 1999a), unbinding forces were found to agree well with the expected logarithmic velocity dependency (Bell, 1978; Grubmüller et al., 1996; Evans and Ritchie, 1997; Isralewitz et al., 1997).

To study the nature and extent of structural heterogeneity of the unbinding pathway, an additional 20 FPMD simulations were carried out and analyzed, each of 300 ps length. These otherwise identical simulations differed in their initial conformations, which were randomly chosen from the last 820 ps of the equilibration phase.

Force profile

View larger version (17K): [in this window] [in a new window]   FIGURE 9   "Force profile" of enforced AN02/DNP-hapten unbinding. Shown is the pulling force exerted on the hapten molecule during unbinding as a function of cantilever position zcant. The force maximum at zcant = 4.0 Å is interpreted as the unbinding force.

Initial unbinding steps

View larger version (23K): [in this window] [in a new window]   FIGURE 10   Three snapshots (a-c) of the initial steps (part A) of the hapten (gray) unbinding pathway. Shown are relevant residues of the binding pocket and the main interactions observed in the simulations such as hydrogen bond and van der Waals contacts (dashed lines). The respective cantilever positions zcant are indicated at the right. The black arrow indicates the pulling direction.

Structural heterogeneity

View larger version (30K): [in this window] [in a new window]   FIGURE 11   Snapshots taken from FPMD simulations with varying unbinding pathways. The left snapshots illustrate the "light chain route," the right ones the "heavy chain route." The cantilever positions of the snapshots are zcant  0 Å (top), zcant  5 Å (middle), zcant  10 Å (bottom), respectively.

View larger version (27K): [in this window] [in a new window]   FIGURE 12   Interaction "fingerprints" (top) and force profiles (bottom) obtained from two representative simulations that took the "light chain route" (left) and the "heavy chain route" (right), respectively. The "fingerprints" show the strength of individual interactions between the hapten molecule and the relevant binding pocket residues "R" and the total interaction energies between the hapten molecule and the light chain [(lightchain)] and the heavy chain [(lightchain)], respectively, as a function of cantilever position zcant. The thickness of the lines denotes the interaction energy.

View larger version (54K): [in this window] [in a new window]   FIGURE 13   Unbinding pathways for 20 FPMD simulations, each with pulling velocity 5 m/s, but slightly differing initial conditions. Shown are the traces (colored lines) of the hapten center of mass and the starting structure (hapten: gray; binding pocket residues: orange). For the two extreme routes, the unbound hapten molecule as seen at the end of the simulation at zcant = 15.0 Å is shown in yellow and green. The red dashed line indicates the force maximum that separates the unbinding pathway into the two parts A and B referred to in the text.

Entropic contribution to the AN02/DNP-hapten bond

View larger version (22K): [in this window] [in a new window]   FIGURE 14   Entropic contribution TSm(zcant) (at T = 300 K) to the free energy landscape along the unbinding reaction pathway as a function of cantilever position zcant. The individual traces have been obtained by taking the m = 50, 55, 60, ... , 100 most significant degrees of freedom into account.

Enforced unbinding for two AN02 mutants

As already discussed, certain residues of the AN02 binding pocket strongly interact with the hapten molecule (see, e.g., Fig. 12). Particularly strong bonds are formed to the light chain residues, and mutations of those affect the AN02/DNP-hapten binding free energy (Martinez-Yamout and McConnell, 1994). In this section we ask whether such mutations can also alter the unbinding force. To that end, the two critical residues shown in Fig. 15, Tyr-33(L) (green) and Ile-96(L) (red), were chosen as mutation sites, and the respective mutant structures of the AN02 antibody were modeled and subjected to FPMD simulations.

View larger version (68K): [in this window] [in a new window]   FIGURE 15   Location of the two site-directed light chain mutants of the AN02/DNP-hapten complex, Y33F(L) (green) and I96K(L) (red), described in the text.

In the first mutant, Y33F(L), the strong hydrogen bond between the hydroxyl group of Tyr-33(L) and the amino group (N8,H29) of DNP-hapten was deleted by replaci