INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION AND CONCLUSIONS REFERENCES

Antibodies represent one of the most extensively investigated classes of proteins. A wealth of x-ray crystallographic data provides detailed structural information about these types of molecules (Padlan, 1994, 1996). Although recently even entire antibodies could be thoroughly studied in all their structural parts (Harris et al., 1998a), most of the data refer to the antigen-binding fragment (Fab). The Fab carries the binding site formed by the complementarity determining regions (CDR) and is thus of primary interest as it constitutes the "typical" part of a given antibody. It is built of two polypeptide chains, heavy chain (H) and light chain (L), each of which is folded into two distinct immunoglobulin (Ig) domains, the N-terminal variable domain (V_H and V_L, respectively), and the C-terminal constant domain (C_H1 and C_L, respectively).

Although the general structural characteristics of Fab molecules and antibodies are well known, comparatively little information is available about their flexibility and dynamics (see, for example, Nezlin (1990); Brekke et al. (1995)). This is in part due to difficulties with obtaining such information experimentally. The structurally most detailed picture is provided by x-ray crystallography, where indications about large-scale molecular flexibility arise primarily from the comparison of structures solved under different conditions (cf. Davies and Chacko (1993); Wilson and Stanfield (1994)), e.g., Fabs crystallized in two or more forms (cf. Sheriff et al. (1987)), Fabs present in the same asymmetric unit of a single crystal (cf. Prasad et al. (1988); Rini et al. (1993)), or the two Fab arms of an intact IgG molecule (cf. Harris et al. (1997, 1998b)). Of primary interest are comparisons between the unliganded and the complexed state (e.g., Stanfield et al. (1993); Rini et al. (1992)), because the corresponding conclusions may reveal motions of functional relevance. Comparisons among x-ray structures have, however, some limitations as far as dynamics are concerned: first, the crystalline environment may not be sufficiently representative for the situation in solution, and second, a comparison between essentially static images representing average structures cannot be expected to provide detailed dynamical information (Petsko, 1996).

A valuable complementary tool to investigate the dynamical properties of proteins in solution is provided by computational molecular dynamics (MD) simulations (Karplus and McCammon, 1983; Karplus and Petsko, 1990; van Gunsteren and Berendsen, 1990; van Gunsteren et al., 1995). They allow following molecular motions at the atomic level and thus help to reveal molecular characteristics otherwise hardly accessible. The method is well established and has been tested and applied on many different systems, including antibody Fab and Fv fragments (e.g., Tanner et al. (1992); de la Cruz et al. (1994); Lim and Herron (1995)). A significant problem, however, is associated with the time scale of the molecular motions that are accessible by simulation, as system size and computational resources impose certain limits. In this context Fab molecules (which consist of >400 amino acids) are already considered as large proteins, and consequently simulation times in the past never exceeded 200 ps. This obviously limits the relevance of statements about motions covering larger scales in time and space.

In this work we present four different MD simulations of significantly extended length for a special Fab. The system of interest is NC6.8, an antibody directed against the sweet-tasting compound NC174 (N-(p-cyanophenyl)-N'-(diphenylmethyl)-N"-(carboxymethyl)guanidine, cf. Fig. 1 in the Methods section). The Fab of this antibody has been analyzed by x-ray crystallography both with and without the hapten, and comparisons between the two structures showed differences of previously unknown magnitude in the domain orientation, most clearly illustrated by a change in the elbow angle of over 30° (Guddat et al., 1994, 1995). In addition, the changes in the binding site provided evidence for an induced-fit type mechanism of ligand binding. The unique character of the observations led to hypotheses about transmitted conformational changes and intramolecular signaling upon complexation. The conclusions, however, were drawn from a single observation and a comparison between two structures in a crystalline environment, which let them appear disputable (Wilson and Stanfield, 1994; Guddat et al., 1995). As a matter of fact, the dynamics of antibody molecules in solution are still not sufficiently understood, and it is still rather unknown how far the changes can go that a hapten induces in an antibody.

In a preliminary account of this work (Sotriffer et al., 1998) two simulations of Fab NC6.8 have already been presented, both starting from the uncomplexed x-ray structure, where in one case the hapten had been introduced into the binding site as initial perturbation. These simulations are further analyzed here. In addition, two new simulations of 700 ps length are presented, both starting from the complexed x-ray structure of NC6.8, in one case complete with the hapten, in the other with the hapten removed. For all four simulations an overall analysis of the quaternary structural features, as well as of the binding site region, is provided.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION AND CONCLUSIONS REFERENCES

Antibody Fab NC6.8 was investigated by means of MD simulations, using four different initial configurations as starting points. They are derived from the x-ray structures of the free, uncomplexed Fab (PDB ID 1cgr), referred to with code "f," and of the complexed Fab (PDB ID 2cgs), referred to with code "c." Depending on whether the hapten is present (code "p") or absent (code "a") in the binding site, the four different simulation systems are given as "fa," "cp," "fp," and "ca." For the simulations fa and cp the corresponding x-ray structure could be used directly. The two additional starting points were generated by modifying these structures with respect to their state of complexation. On the one hand, the hapten was inserted into the binding site of the free Fab structure (simulation fp); on the other hand, the ligand was removed from the binding site of the complexed Fab structure (simulation ca). The insertion of the ligand for simulation fp was done with the aid of interactive modeling tools (SYBYL 6.3, Tripos Associates (1996)) and guided by the position of the ligand in the x-ray structure of the complex; conformational changes in the antibody itself were not applied. The modifications for the removal of the ligand to obtain simulation system ca consisted simply in deleting the ligand's coordinate entries in the coordinate file of the complex.

All four systems were set up for simulation in an identical way, except for differences arising from the presence or absence of the ligand. All preparing steps and simulations were carried out with the AMBER 4.1 suite of programs (Pearlman et al., 1995a,b). Because the all-atom AMBER force field of Cornell et al. (1995) was used for the simulations, missing hydrogen atoms had to be added to the starting coordinates. This resulted in a total of 6491 atoms for the 433 amino acid residues of the Fab and 49 hapten atoms. Both the Fab and the ligand had a net charge of zero.

The setup of the ligand, hapten NC174, required the explicit assignment of suitable atom types, which are shown in Fig. 1. Additional force field parameters required for structural units not covered by the standard AMBER force field are listed in Table 1. They were derived by comparison with existing parameters, experimental data, and results from ab initio calculations. Following the philosophy of the Cornell et al. (1995) AMBER force field, charges for the ligand were calculated by fitting to the HF/6-31G* electrostatic potential. The corresponding ab initio calculations were performed with GAUSSIAN94 (Frisch et al., 1995), the restrained electrostatic potential fit with the RESP program (Bayly et al., 1993; Cornell et al., 1993).

View larger version (8K): [in this window] [in a new window]   FIGURE 1   Atom types used for NC174. Hydrogen atoms of the aromatic rings are not shown explicitly; they were all assigned atom type HA.

To relax internal strains and remove bad initial contacts, the starting structure of each system (fa, fp, cp, ca) was subjected to a short energy minimization consisting of 20 steps steepest descent and 280 steps conjugate gradient, using an effective distance-dependent dielectric constant of  = r. Subsequently, the Fab or Fab-ligand complex was placed in a rectangular box of TIP3P water molecules (Jorgensen et al., 1983), with a minimum solute to wall distance of 8 Å. This resulted in average box sizes of 96 Å × 74 Å × 62 Å and ~10,500 water molecules in each system. The solvated systems were again energy-minimized (50 steps steepest descent, 450 steps conjugate gradient;  = 1), primarily to optimize the interactions at the protein-water interface.

The simulations were started from the minimized and solvated systems. For equilibration, the protein (and ligand) atoms were initially frozen in their position and only the water molecules were allowed to move. Under this condition, the system was heated to 300 K over 4 ps and subsequently cooled down to 80 K in 1 ps. Then the protein atoms were allowed to move as well and the system was heated to 300 K over 15 ps of simulation time. From this point on, the simulation was carried on at 300 K and 1 bar (NPT conditions). The temperature was kept constant by coupling to a heat bath through the Berendsen algorithm (Berendsen et al., 1984) using separate solute and solvent scaling. Pressure was adjusted by isotropic position scaling using a Berendsen-like algorithm. Covalent bonds to hydrogen atoms were constrained by the SHAKE algorithm (Ryckaert et al., 1977). A time step of 2 fs was used. The simulations were carried out under periodic boundary conditions, and a twin-range residue-based cutoff of 8 and 10 Å was applied to the non-covalent interactions; the pairlist was updated every 10 time steps (0.02 ps). For analysis, energy data were saved every 0.02 ps, coordinates every 0.1 ps. Each simulation presented here was carried out for 700 ps.

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION AND CONCLUSIONS REFERENCES

Simulation fa: free Fab structure, hapten absent

Energetical equilibration in the simulation of the uncomplexed Fab required ~100 ps. In general the analyses and trajectory averages refer therefore to the period from 100 to 700 ps. The simulation produced a sufficiently stable trajectory, representative of an equilibrium state in solution. The root-mean-square (rms) deviation of the C atoms from the starting x-ray structure amounts to 3.09 (±0.11) Å for the last 100 ps, which is comparable to other Fab simulations (Lim and Herron, 1995). If only the C atoms of the structurally most conserved central core of each Ig domain are used for superposition and rmsd calculation (40 residues in each variable domain, 35 residues in each constant domain), a value of 2.23 (±0.15) Å is obtained, which shows the same progression as the all-C rmsd. Both values suggest interdomain adjustments and movements, because for the single domains significantly lower values (~1 Å) are measured.

Quaternary structural dynamics were first analyzed with respect to the so-called Cys-Trp-Cys triads. These are highly conserved structural elements present in each Ig domain and may be used to define central points and planes as reference for distance and angular measurements to evaluate the relative domain orientation. As shown in detail in Sotriffer et al., 1998, the trajectory averages reproduce the experimental values determined for the uncomplexed x-ray structure very well. The rms fluctuations for the distance averages are between 0.4 and 0.8 Å, while the angular measurements show fluctuations of 5.6-6.8°. These values indicate moderate movements of the domains with respect to each other.

Parameters often used to characterize the domain orientation in Fab molecules are the pseudodyad angles and the elbow angle. The former is defined by a rotation of the V_L (or C_L) domain around a pseudodyad axis to achieve optimal superposition with the V_H (or C_H1) domain (the C atoms of the conserved cores were used for superposition; the calculations were done with routines of the ALIGN program (Satow et al., 1986)). The angle formed by the V pseudodyad axis and the C pseudodyad axis is the elbow angle. The simulations yield average values for the pseudodyad angles that are slightly lower than in the x-ray structure, with rms fluctuations of 2-3°. More importantly, the elbow angle (experimental value 189.3°) is well reproduced by the trajectory average of 187.8° and its rms fluctuations of 2.0° (minimum value 183.4°, maximum value 194.2°). Most interesting is the time course of the elbow angle, as it suggests a periodic hinge-bending motion. As revealed by the corresponding autocorrelation function and its Fourier transform, this motion appears to be a vibration with a period of 137 ps, onto which a faster fluctuation of smaller amplitude and 23-ps period is superimposed.

This elbow motion suggests time-correlated protein domain motions of the variable domains with respect to the constant domains. To investigate this further, a cross-correlation analysis of atomic displacements was performed to obtain so-called dynamic cross-correlation matrices (DCCM) (Ichiye and Karplus, 1991; Swaminathan et al., 1991). These are calculated by the expression

(1)

where r_i is the displacement from the mean position of the ith amino acid (the coordinate of an individual amino acid was calculated from the mean of the N, C, and C backbone atom coordinates). In a contour plot of the matrix [C_ij] strong correlations in atomic motion show up as large off-diagonal crosspeaks. As pointed out by Arnold and Ornstein (1997), attention has to be paid to the reference frame used for structural superposition of the MD snapshot structures. Typically, a global superpositioning using all backbone atoms is used. The DCCM shown in Fig. 2 was obtained this way. However, as described in detail by Arnold and Ornstein (1997), the use of a global reference frame can potentially obscure time-correlated domain motions of the protein. Therefore, two further DCCM analyses were carried out, one using the backbone of the variable domains as reference frame (Fig. 3), the other using the backbone of the constant domains (Fig. 4). The effect can be seen immediately. While for the global reference frame high correlations are observed only within the domain (mostly along the characteristic -sheet structure), the separate reference frames clearly reveal highly correlated motions of the entire domain pair. If the variable domains serve as stationary points of reference, the C_H1-C_L domain pair emerges as a concertedly moving unit. Conversely, if the constant domains are used as reference, the Fv displays high correlations. These results support the view of the elbow motion as a concerted vibration of the variable domains with respect to the constant domains. Also, the time interval of 200 ps used for the DCCMs presented here fully covers the period of the elbow vibration as suggested by the autocorrelation analysis.

View larger version (53K): [in this window] [in a new window]   FIGURE 2   Dynamic cross-correlation map of the time-correlated backbone atom motions of trajectory fa, calculated for the 200-400-ps time interval, using the entire structure (all 433 residues) as coordinate reference frame. Only cross-correlations >0.25 are shown. Positive correlations are shown in the lower triangle, negative correlations in the upper triangle. The residues are numbered sequentially, starting with the H chain (1-214) and continuing with the L chain (215-433).

View larger version (59K): [in this window] [in a new window]   FIGURE 3   Dynamic cross-correlation map of the time-correlated backbone atom motions of trajectory fa, calculated for the 200-400-ps time interval, using the variable domains (Fv, residues 1-113 and 215-322) as coordinate reference frame. Other explanations as for Fig. 2.

View larger version (53K): [in this window] [in a new window]   FIGURE 4   Dynamic cross-correlation map of the time-correlated backbone atom motions of trajectory fa, calculated for the 200-400-ps time interval, using the constant domains (CH1 and CL, residues 121-214 and 330-433) as coordinate reference frame. Other explanations as for Fig. 2.

Special rms and distance measurements combined with visual inspection served to analyze the structural and dynamical features of the CDR region and the immediate binding site. The average rms deviations for the CDRs are shown in Table 2. The values in column A are based again on the superposition of all C atoms of the Fab. For comparison: the corresponding average (100-700 ps) for the entire C structure is 2.63 (±0.32) Å. Therefore, CDR H2, L1, and L3 show larger deviations than the overall structure, while CDR H1, H3, and L2 show smaller ones. The fluctuations, however, are larger for all six CDRs. If only the residues of the CDRs are used for superposition (column B in Table 2), somewhat smaller values are obtained, as expected. Interestingly, the superposition of the single CDRs (column C) yields significantly smaller values (all below 1 Å). The backbone conformations of the single CDR loops thus appear to be largely conserved, while their relative orientation is subject to alterations.

View this table: [in this window] [in a new window]   TABLE 2   Simulation fa average C rms deviations from the x-ray structure for the single CDRs, using different sets of atoms as reference for superposition and evaluation: all C atoms of the Fab (column A), the C atoms of all CDRs (column B), and the C atoms of each single CDR loop (column C). All values are given in angstroms

Of interest is a comparison with C rms values of the CDRs provided by de la Cruz et al. (1994) for another antibody. In considerably shorter simulation times (75 ps after 50 ps equilibration) they find for the single CDR loops of the uncomplexed antibody rms deviations between 0.54 and 1.45 Å, and only CDR L3 (0.54 Å) and CDR H1 (0.61 Å) show values below 1 Å. Thus, the CDR rms deviations are higher than for NC6.8, although the order of magnitude suggests a similar intrinsic variability of the CDR loops.

As far as the immediate binding site is concerned, the analysis is focused on four aromatic amino-acids that flank the binding site and somehow form its corners: Tyr-L32 (L1), Tyr-L96 (L3), Trp-H33 (H1), and Tyr-H96 (H3) (residues numbered according to Kabat et al. (1991); cf. also Fig. 7). To judge their mobility and the size-variability of the binding site, distance measurements were carried out, both for the backbone (C atom distances) and for the side chains (distance between the geometric centers of the aromatic rings). The results are shown in Table 3.

The average values generally remain close to the distances measured in the x-ray structure. However, considerable fluctuations are observed for the side chains, which is most likely a consequence of their exposed position. As shown also by visual inspection, Tyr-L32 and Tyr-L96 remain essentially unchanged. In contrast, the indole side chain of Trp-H33 is slightly reoriented. The Tyr-H96 side chain remains in the "down" orientation (cf. Guddat et al. (1994)) and does not move upward, as experimentally observed in the complex. However, the H3 loop moves a little toward the binding site center and appears to reduce the size of the entrance thereby. Accordingly, nearly all side chain distances are on average smaller than in the experimental structure.

In summary, the results for simulation fa appear reasonable and suggest that the dynamics of the Fab can be described by the applied simulation method. The average values of many parameters reproduce the experimental data for the crystal structure fairly well, which indicates a sufficiently stable simulation. Time courses and fluctuations of the structural parameters can therefore be expected to provide relevant insights about the flexibility of this antibody fragment.

Simulation fp: free Fab structure, hapten present

How does the simulation of the uncomplexed structure proceed if the hapten is inserted into the binding site? As expected, the insertion of the ligand leads to a perturbation of the system, which is reflected by a significantly prolonged equilibration phase. The potential energy requires ~200 ps to reach an equilibrated value. If not stated otherwise, the following analyses therefore refer to the period from 200 to 700 ps.

The rms deviations from the starting structure (x-ray structure of the uncomplexed Fab) indicate that considerable changes occur: the average rmsd for all C atoms over the last 100 ps is 6.24 (±0.15) Å. In a "standard" simulation, a value of this size would indicate questionable quality and stability. However, in the case of simulation fp, the purpose is not to generate an equilibrium ensemble, which on average corresponds to the starting structure. The objective is rather to follow which structural changes are induced by the perturbation. In this sense, the rms deviations do not serve as quality criteria, but as first indications of significant structural changes. Interestingly, similar rms values are obtained if the x-ray structure of the complex instead of the uncomplexed structure is used as reference. In some cases even lower values result with respect to the experimental structure of the complex. Most significant is this difference for the L chain: using the core residues of the V_L and C_L domains for the C rms fit, a trajectory average (200-700 ps) of 4.16 (±0.33) Å is obtained with respect to the free x-ray structure, and 3.05 (±0.18) Å with respect to the complexed x-ray structure. The corresponding values for the last 100 ps are 4.23 (±0.17) Å and 3.13 (±0.12) Å, respectively. This is a first indication that the Fab with the inserted hapten proceeds toward a state that resembles the x-ray complex structure at least as much as the starting uncomplexed structure.

A more detailed structural analysis reveals that indeed significant quaternary structural changes take place which in many aspects correspond to the experimental differences between free and complexed x-ray structure. The elbow angle rapidly decreases during equilibration and reaches an average value of 168.9° for the last 500 ps, with a maximum of 176.0° and a minimum of 160.7° (the experimental value for the complex is 152.4°, as compared to 189.3° for the uncomplexed Fab). The L-chain is elongated, while the H-chain becomes more flexed. This is seen both in the Cys-Trp-Cys triad distances between the V_L and C_L domain and between the V_H and C_H1 domain, as well as in the end-to-end distances measured between the C atoms of the N-terminal and the C-terminal residues of the H-chain and the L-chain, respectively (see Sotriffer et al. (1998) for detailed values).

As far as the triad distances between the two variable domains and between the two constant domains are concerned, experimentally there is no real difference between the complexed and the uncomplexed state. In the fp simulation, however, these distances are somewhat increased and trajectory averages of ~25 Å are obtained instead of the experimental ~21 Å. The ligand thus seems to penetrate into the interface between the V_H and the V_L domain, a process that has been mentioned as the possible initiating event for the structural changes (Guddat et al., 1995).

To analyze the situation in the binding site around the hapten and compare it with the fa simulation, distance measurements between the four essential aromatic residues were again carried out (cf. Table 4). The average distances between the C atoms show in part better correspondence with the complexed crystal structure than with the uncomplexed one. However, the difference between the complexed and uncomplexed state is rather small (~0.5 Å) and in the same order of magnitude as the fluctuations of the simulation values. Therefore, the averages do not serve to illustrate a transition between complexed and uncomplexed state, but rather to show that key features of the binding site architecture remain intact after insertion of the hapten. The fluctuations of the distances are (with one exception) smaller or equal in size as in simulation fa. This is also observed for many other parameters not further discussed here. It reflects the expectation that the flexibility of the binding site is slightly reduced after ligand binding and a more compact unit is formed together with the hapten.

Also, the side chains show in sum smaller fluctuations than in the uncomplexed simulation. However, the averages do in general not compare very well with the experimental data for the complex. As shown by a visual analysis, the side chains of Tyr-L96 and Trp-H33 are oriented as observed in the complex (Tyr-L96 in close contact and parallel to the cyanophenyl group of NC174, Trp-H33 in interaction with the central guanidinium unit). In contrast, the side chain of Tyr-L32 moves slightly away from the binding site. Most importantly, the swing (upward movement) of Tyr-H96 does not occur. Instead, this side chain remains in its initial position and points away from the binding site without showing significant contacts with the ligand.

The ligand itself appears generally well accommodated within the binding site. However, a significant difference to the experimental complex structure can be observed at one side of the binding pocket. While normally the contacts between CDR L3 and CDR H2 form the border of the binding site, the distance between these two CDRs somewhat increases in the simulation, and an opening of the binding site arises. While this has no direct consequences for the position of the ligand, some interactions with CDR H2 are obviously lost. Apparently the strains after the initial insertion of the ligand are sufficient to lead to an opening on this part of the binding site.

Good agreement is found for the shortening of the distances between Gly-L91 and CDR H3 and between Tyr-H96 and Ser-H98. The simulation gives average distances of 7.7 (±0.5) Å between O(Gly-L91) and C(Tyr-H96), 4.5 (±0.4) Å between O(Gly-L91) and C(Ser-H97), and 6.7 (±0.5) Å between O(Gly-L91) and C(Ser-H98). The corresponding transitions from the uncomplexed to the complexed x-ray structure are in the first case 8.6  8.4 Å, in the second 7.2  5.5 Å, and in the third 8.7  7.2 Å. As in the experiment the opening of the -turn between Tyr-H96 and Ser-H98 is observed, with the loss of the hydrogen bond between the amide of Tyr-H96 and the carbonyl of Ser-H98. The substitution of this hydrogen bond by a hydrogen bond between Tyr-H96 (carbonyl-O) and the guanidinium group of NC174 is not observed. This is also due to the fact that the upward movement of the Tyr-H96 side chain does not occur and no further contact with the hapten is achieved. Other aspects of the interaction with NC174 agree well with the experimental findings, as for example the distances to Tyr-L96 and Trp-H33. Furthermore, the hapten almost exactly retains its conformation during the simulation, with trajectory averages for the torsion angles closely matching the experimental starting values.

In summary, not all details of the complexation and its structural consequences are reproduced as in the experimental x-ray structures, but nevertheless the agreement in many aspects is striking. It is in fact remarkable that in a simulation of 700 ps the simple insertion of the hapten into the binding site causes the Fab to achieve a significantly altered quaternary structure (not only local changes), which resembles more the experimental complexed state than the uncomplexed x-ray structure.

Simulation cp: complexed Fab structure, hapten present

For comparative purposes, a simulation using the experimental complex structure as starting point was carried out as well. It was expected to observe a similar behavior as for simulation fa, i.e., fluctuations around an equilibrium state which, on average, corresponds fairly well to the x-ray structure. The energetical equilibration required again ~100 ps, with the total potential energy showing no further drift and fluctuations of ~0.1% after 100 ps.

Given the findings in simulation fa and the energetical behavior of simulation cp, the rms measurements shown in Fig. 5 provided a big surprise. While for the first 100 ps after energetical equilibration the all-C rms deviation from the starting x-ray structure is on average 2.3 Å, the value increases considerably during the following 500 ps, reaching ~5.5 Å in the last 100 ps. These deviations result from changes in the overall structure because a similar curve with slightly lower values is obtained when only the conserved core residues are used (rms2 in Fig. 5). The tertiary structure of the single domains is not affected, as shown by the superpositions and rms deviations of the core C atoms of the single domains. For V_L, for example, the trajectory average (100-700 ps) is 0.71 (±0.07) Å, which is even lower than the 0.78 (±0.06) Å in the case of simulation fa. Similar values are obtained for the other domains (V_H, 0.83 ± 0.16 Å; C_L, 0.71 ± 0.07 Å; C_H1, 0.77 ± 0.10 Å), whereby the values for the H-chain domains are slightly larger than for the L-chain domains. First indications about the changes occurring in the quaternary structure were obtained by calculating rms deviations for the H- and the L-chain separately, which showed that the deviations for the H-chain are much larger (by ~3 Å) than for the L-chain.

View larger version (22K): [in this window] [in a new window]   FIGURE 5   Rms deviations from the x-ray structure of the complex, measured in simulation cp for all C atoms (rms1), for all C atoms of the central structurally conserved domain cores (rms2), and for the central core C atoms of the VL domain (rms3).

This observation was confirmed by a series of other analyses. Distance measurements between the C atoms of the terminal residues of both chains show that during the simulation the end-to-end distance for the H-chain is significantly reduced, while the values for the L-chain fluctuate around an average somewhat below the experimental value. Starting at the initial 48.1 Å, the distance between the H-chain termini decreases to an average of 31.8 (±1.3) Å for the last 100 ps. In contrast, the L-chain end-to-end distance shows a trajectory average of 61.7 (±1.0) Å (last 100 ps: 61.9 ± 0.8 Å), which corresponds to a slight shortening compared to the experimental 67.0 Å. It thus appears that for the H-chain the trend toward shorter distances and a closure of the V_H-C_H1 interface (observed experimentally for the transition from the uncomplexed to the complexed state) is persisting.

The measurements carried out using the Cys-Trp-Cys triads confirm that the changes are essentially restricted to the H-chain and to the orientation of the V_H domain relative to the C_H1 domain. The triad distance between these two domains decreases by ~10 Å during the simulation and reaches a rather stable value of 24.4 Å in the last 100 ps. The other average distances correspond fairly well to the experimental values and show no change during the simulation. A similar situation is seen for the angles, where the V_H-C_H1 angle changes from 49.6° to 102.3° (the fact that a larger angle results, although the movement of the domains is described as "closure," is due to the orientation of the triad planes). However, for the other angles deviations are observed as well, but only in the range of 10-20°. The angles react more sensitively to structural changes, and obviously a drastic change in the H-chain cannot occur without certain adaptations in the orientation of the other domains (especially V_L-C_L).

In accordance with the observations described so far, rather drastic changes are also observed for the elbow angle (cf. Fig. 6). Starting at the experimental 152.4°, the angle decreases to 135° during equilibration and finally reaches a value of 106°, which is maintained on average during the last 100 ps. According to a visual inspection of the quaternary structure, it appears that an elbow angle of this size corresponds to the physical limit of the Fab bend (to our knowledge, the smallest experimentally observed value is 127° for antibody 8F5 (Tormo et al., 1994)). The reason is that already new van der Waals contacts between the variable and the constant domains appear, which seem to hinder a further bending. Indeed, in the last 100 ps twice as many van der Waals contacts between V_H and C_H1 are observed as in the initial stages of the simulation.

View larger version (41K): [in this window] [in a new window]   FIGURE 6   Backbone drawings of Fab NC6.8 complexed with hapten NC174: the experimental starting structure (left) and the trajectory snapshot after 700 ps of simulation cp (right).

As far as the binding site is concerned, the question is whether similarly drastic local changes occur. Because distances and angles between the V_H and V_L domain change only comparatively little, larger changes within the CDRs must not necessarily be expected. To investigate this and allow for comparisons with the simulations fa and fp, rms deviations for the CDR loops and distances between binding site residues were again calculated. As far as the rms deviations are concerned, a fit with respect to all C atoms of the Fab was not done, because it does not appear reasonable in this case. Therefore, Table 6 shows only the results of the superposition of all CDR C atoms (column B), and of the C atoms of each single CDR loop (column C).

View this table: [in this window] [in a new window]   TABLE 6   Simulation cp average C rms deviations from the x-ray structure for the single CDRs, using different atom sets as reference for fitting and evaluation: the C atoms of all CDRs (column B), and the C atoms of each single CDR loop (column C). All values are given in angstroms

The values confirm that no drastic changes occur. Of special interest is a comparison with the results obtained for simulation fa. It shows that except for CDR H1, the rms deviations in simulation cp are smaller than in simulation fa. In most of the cases also the fluctuations are smaller than in the uncomplexed state. The ligand thus seems to stabilize the conformations of the CDR loops, a finding that corresponds to the expectations and could also be seen in other antibody simulations (de la Cruz et al., 1994).

The slightly larger deviations and fluctuations for the CDR loops of the H-chain are most probably related to the quaternary structural rearrangements of the H domains. This is further supported by the fact that these rms curves show a similar trend as other parameters that characterize the overall structure: an increase during the central simulation phase and a stabilization at higher levels toward the end. The visual analysis indicates that due to the differences in the V_H-V_L orientation (triad distance 1.5 Å larger, triad angle 16° larger) CDR H1 and H2 move somewhat away from the ligand and the binding site center. Because CDR H3 remains in close contact with the ligand, however, a larger distance between CDR H3 and CDR H1/H2 results.

View larger version (38K): [in this window] [in a new window]   FIGURE 7   The binding site of Fab NC6.8 with hapten NC174 displayed in dark gray. The figure shows the experimental conformations that served as starting point for simulation cp and simulation ca. The image was generated with MOLSCRIPT (Kraulis, 1991).

These observations are confirmed by the binding site distance measurements shown in Table 7. All the distances to Trp-H33 are increased (with exception of the side chain distance Trp-H33-Tyr-H96, where the side chain orientation compensates the increase in the backbone distance). Also, the fluctuations are larger than for the other distances. This is due to the fact that these parameters show again the typical time course of other H-chain parameters (increasing trend in the central simulation phase).

The other values show a better accordance with the experimental references and the fluctuations are generally even smaller than the corresponding values in simulation fa. Other structural units of the binding site are relatively well conserved, too, such as the distances between Gly-H91 and CDR H3 (Tyr-H96, Ser-H97, Ser-H98), which increase only by some tenths of an angstrom. As far as the position of the ligand within the binding site is concerned, its tight fit to CDR L1, L3, and H3 is well conserved, while the distance to CDR H1 (Trp-H33) is increased. The contacts with Tyr-L96 and Tyr-H96 remain intact, and the average distances correspond almost exactly to the experimental values (Tyr-H96 itself remains in the upward oriented position, as in the starting structure). Also, the hydrogen bond between Tyr-H96 and the guanidinium of NC174 is conserved most of the time. As far as the conformation of the ligand is concerned, most of the torsion angles remain on average close to their starting value. Somewhat different torsional preferences are observed only for the diphenylmethyl unit, which is not surprising given its exposed nature.

An obvious question concerning simulation cp is of course how to interpret the results with respect to the experimental facts from x-ray crystallography and the simulation methodology. This issue will further be dealt with in the Discussion section.

Simulation ca: complexed Fab structure, hapten absent

As the fourth variant, simulation ca served the purpose to investigate whether the removal of the hapten from the binding site would cause the reverse change as its insertion. Because in structural terms the applied perturbation is definitely smaller (compared to simulation fp), the probability to observe such changes must expected to be lower. As a matter of fact, the energetical equilibration required about the same time as in simulation fa and cp, but not as much as in simulation fp.

The rms deviations from the x-ray structure of the complex reveal a rather surprising behavior, given the observations made so far: Fig. 8 shows that the quaternary structure remains stable along the trajectory. The trajectory averages (100-700 ps) are even below the values measured for simulation fa. The average rms deviation is 2.52 (±0.21) Å for all C atoms and 1.91 (±0.22) Å if using only the C atoms of the "core residues." Also for the single domains very low values are obtained: V_L 0.69 (±0.08) Å, V_H 0.60 (±0.06) Å, C_L 0.74 (±0.08) Å, and C_H1 0.82 (±0.14) Å.

View larger version (21K): [in this window] [in a new window]   FIGURE 8   Rms deviations from the x-ray structure of the complex, measured in simulation ca for all C atoms (rms1), for all C atoms of the central structurally conserved domain cores (rms2), and for the central core C atoms of the VL domain (rms3).

The fact that the initial quaternary structure is generally well maintained is also reflected by the distance measurements between the Cys-Trp-Cys triads, which give average values that reproduce the experimental values for the complex. The angular measurements, however, suggest that nevertheless after removal of the hapten some small adjustments occur in the domain orientations. This is also shown by the elbow angle, which is on average 140° and fluctuates between 135° and 145°. Compared to the starting structure it is thus reduced by ~10°, which means that it shows an adaptation to a lower value instead of a higher one, which would correspond to the experimental uncomplexed state.

As far as the CDRs are concerned, no important changes are observed as well. The rms deviations for the C atoms of the six CDR loops suggest that in simulation ca these regions are at least as well conserved as in simulation cp where the hapten is present in the binding site. The fluctuations are of the same size and the deviations after superposition of the entire CDR region are even somewhat smaller (with exception of CDR L3). In any case, no real perturbation can be observed after removal of the hapten.

Regarding the four aromatic amino-acids of the binding site, however, the removal of the hapten has some consequences for the position of the side chains. As revealed by distance measurements and shown by visual inspection, the Tyr-L96 side chain (in the complex oriented parallel to the cyanophenyl ring of NC174) tilts into the binding site, while Trp-H33 turns away toward the top. This way the distances between Tyr-L96 and Tyr-L32/Tyr-H96 become smaller, those to Trp-H33 somewhat larger than in the complexed state. The Tyr-H96 side chain remains in the "upward" directed conformation, where it is completely solvent-exposed and has considerable conformational freedom in absence of the hapten. Accordingly, large fluctuations are observed for the distances to the Tyr-H96 side chain.

Apparently the simple removal of the ligand does not produce a sufficient perturbation that would result in significant structural changes or even transitions toward the experimental uncomplexed state within the simulation time. Interestingly, not even those changes are observed that occur in simulation cp in the presence of the ligand. Rather, the starting structure of the complex is largely maintained, with lower rmsd values than in the "unperturbed" simulation cp.

	DISCUSSION AND CONCLUSIONS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION AND CONCLUSIONS REFERENCES

The simulations presented here belong to the most extensive simulation studies on antibody systems hitherto reported, with respect both to the simulation length and the comparison between four different starting conditions. This enables us to obtain better estimates of the flexibility of Fab molecules in solution and to complement experimental findings derived from comparisons between "static" structures. In addition, the comparison of the four different simulations also allows us to judge the quality, reliability, and limits of the applied methodology.

Simulation fa (starting from the uncomplexed Fab NC6.8) may be regarded as a "standard" simulation and can be used to evaluate the quality of the technique. As revealed by comparisons with other simulation studies of Fab (Lim and Herron, 1995) and Fv molecules (de la Cruz et al., 1994) and according to common standard methodology, simulation fa completely fulfills the general quality criteria. Although the other cited simulations were performed in a similar way (explicit solvent, periodic boundary conditions, cutoff for non-covalent interactions) and for significantly shorter simulation times (below 200 ps), on average lower rms deviations from the x-ray structure were obtained in simulation fa. Eventual limitations should therefore not result from the special procedure used here, but rather be related to the general problems of MD simulations.

Based on the observation that on average the x-ray structure is well maintained in simulation fa and that essential structural parameters are well reproduced by the trajectory, statements about the flexibility of this Fab should be sufficiently reliable. Of special interest in this context is the elbow angle, the most important parameter used to characterize the Fab quaternary structure. In general, the experimentally observed elbow angle of Fab molecules in different states varies by <15°. The variation in the elbow angle over a range of 11° observed in simulation fa suggests that the differences known from experimental data need not to be caused by drastic forces, but may rather correspond to different states which, due to the flexibility of the molecule, may also be occupied in equilibrium in solution. Of interest is the observation of a vibration-like behavior of the variable domains with respect to the constant domains and the fact that this is observable in a simulation of 700 ps. A period of ~140 ps suggests that in fact the domains are moving comparatively fast. A fast relative movement can be reasonable because only very few contacts between the two domain pairs exist, allowing them to move rather unhindered.

Simulation fp was started from a perturbed initial state where the ligand was inserted into the binding site of the unperturbed Fab structure. The resulting trajectory is remarkable: the simple insertion of the small ligand into the binding site causes significant changes in the quaternary structure, within comparatively short time scales. The observed movements clearly exceed those observed for the equilibrated state of simulation fa. Furthermore, it is interesting that the changes of many structural parameters tend toward the experimental values measured for the complexed state. Of course this does not apply to all structural aspects, but the expectation of an exact reproduction of the experimental complex would be rather unrealistic. Within the binding site, for example, the upward flip of Tyr-H96 does not occur, which may be due to the limited time scale or to the special nature of the starting structure, which may not correctly reproduce the important initiating moment. Even more remarkable is therefore the fact that changes in the quaternary structure correspond to experimental differences between the uncomplexed and the complexed state. The initiating event may be the entering of the ligand into the space between the V_H and V_L domain, as reflected also by the increased distance between these domains in the simulation compared to experiment. In any case the simulation suggests that the ligand may indeed be able to elicit allosteric-like effects. The unusual size of the experimentally observed effects should therefore not primarily be dictated by the crystalline environment, but rather apply to the pure solvated state as well.

In contrast to the two simulations discussed so far, simulation cp reveals a rather unexpected behavior. Starting from the experimental complex structure, a trajectory should be obtained that shows characteristics similar to simulation fa. However, drastic changes result for the quaternary structure (but not for the binding site): the elbow angle is further decreased, the H-chain is further flexed. This corresponds to a continuation of the trends observed experimentally for the transition from the uncomplexed to the complexed state. Apparently, the effect of the ligand on the Fab structure persists even if the simulation is started at the x-ray structure of the complex.

That the characteristics of simulation cp are indeed due to the presence of the ligand is documented by simulation ca. In this case the ligand was removed from the binding site and the complex structure was subjected to simulation without hapten. Surprisingly, the starting structure was fairly well maintained and apart from some smaller changes a structurally stable trajectory was obtained, comparable to simulation fa. The fact that no drastic changes could be observed in the binding site (such as the initial collapse in the study of Lim and Herron (1995)) is probably due to the different solvation procedure: here, it was applied after ligand removal, which allowed water molecules to fill the space from the beginning of the simulation, with the result that a larger perturbation was avoided.

The results of the four simulations may therefore be summarized as follows: the presence of the hapten induces structural changes in the quaternary structure, while in the absence of the ligand fluctuations around an average equilibrium state largely corresponding to the starting structure are obtained. Thus, the hapten shows clear effects and the Fab structure in solution reacts very sensitively to its presence. An important point may be that the ligand like a wedge deeply enters the cleft between V_H and V_L, thereby inducing structural changes.

The central question, however, is whether the simulations cp and ca reflect real features. The results of simulation ca still appear rather plausible, considering that no external forces and no large perturbations are present. The Fab thus relaxes toward the nearest stable state and does not show any tendency for larger transitions within the simulation time. Whether this state has a real significance or not is hard to say. Assuming reversibility of the ligand-induced effects, the system should return to its original state in the absence of the ligand. However, experimental details about the reversibility and time scales of the observed transitions in the NC6.8-NC174 system are not known. In contrast, it is well known that proteins are generally characterized by highly complex energy surfaces and may occupy numerous stable or metastable conformational substates (Frauenfelder et al., 1991). In this sense it could be possible that after removal of the hapten Fab NC6.8 for some time (according to the simulation at least some hundred picoseconds) adopts states similar to those found in simulation ca.

An interpretation of simulation cp appears more difficult. Although the simulation in a certain sense confirms the structural effect of ligand binding observed also in simulation fp, it would be expected that starting from an "unperturbed" x-ray structure a trajectory should result that shows only fluctuations around the starting structure. Here, in contrast, the movements caused by the complexation are continued. The only possible interpretation for this observation is that the experimental complex structure is stabilized by the crystal, while in solution a stronger bending is achieved. The plausibility of this assumption is hard to assess, as too little information about Fab molecules in solution is still available. Clearly, the magnitude of the changes leading to an elbow angle of 106° appears rather drastic, considering that the smallest elbow angle currently known from crystal structures is 127°.

In conclusion, one should therefore also think about possible limitations of MD simulations that may lead to artifacts, primarily the ultimately unsatisfying treatment of the electrostatic interactions by a cutoff method. Although the here applied variant of a residue-based twin-range cutoff helps to somehow reduce the problems, it is well known that cutoff methods can lead to artifacts and significantly influence the simulation (Schreiber and Steinhauser, 1992; Saito, 1994). With PME (Particle Mesh Ewald; Cheatham III et al., 1995) an alternative and potentially superior method would be available, which is more and more becoming the standard in MD simulations. However, for systems of Fab size the method is still extremely demanding in terms of currently available computing resources. Furthermore, the periodicity imposed to solutions that are inherently non-periodic systems can give rise to artifacts as well (Hünenberger and McCammon, 1999). For comparative purposes PME simulations for Fab NC6.8 would nevertheless be highly interesting and may hopefully become feasible in the not too distant future.

For the moment no definitive interpretation can be given for the unexpected observations made in simulation cp. It may, however, be argued that with the same simulation protocol reasonable results were obtained for the other three simulations. This could also suggest that simulation cp is not simply an artifact, but rather an indication that the Fab flexibility may go beyond currently known limits.