This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Zweckstetter, M. Articles by Bax, A. Search for Related Content PubMed PubMed Citation Articles by Zweckstetter, M. Articles by Bax, A.

Prediction of Charge-Induced Molecular Alignment of Biomolecules Dissolved in Dilute Liquid-Crystalline Phases

Markus Zweckstetter ^*, Gerhard Hummer  and Ad Bax 

^* Max Planck Institute for Biophysical Chemistry, Am Fassberg, Göttingen, Germany; and Laboratory of Chemical Physics, National Institute of Diabetes and Digestive and Kidney Diseases, National Institutes of Health, Bethesda, Maryland USA

Correspondence: Address reprint requests to Dr. Markus Zweckstetter, Max Planck Institute for Biophysical Chemistry, Am Fassberg 11, D-37077 Göttingen, Germany. Tel.: 49-551-201-2220; E-mail: mzwecks{at}gwdg.de.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
Alignment of macromolecules in nearly neutral aqueous lyotropic liquid-crystalline media such as bicelles, commonly used in macromolecular NMR studies, can be predicted accurately by a steric obstruction model (Zweckstetter and Bax, 2000). A simple extension of this model is described that results in improved predictions for both the alignment orientation and magnitude of protein and DNA solutes in charged nematic media, such as the widely used medium of filamentous phage Pf1. The extended model approximates the electrostatic interaction between a solute and an ordered phage particle as that between the solute's surface charges and the electric field of the phage. The model is evaluated for four different proteins and a DNA oligomer. Results indicate that alignment in charged nematic media is a function not only of the solute's shape, but also of its electric multipole moments of net charge, dipole, and quadrupole. The relative importance of these terms varies greatly from one macromolecule to another, and evaluation of the experimental data indicates that these terms scale differently with ionic strength. For several of the proteins, the calculated alignment is sensitive to the precise position of the charged groups on the protein surface. This suggests that NMR alignment measurements can potentially be used to probe protein electrostatics. Inclusion of electrostatic interactions in addition to steric effects makes the extended model applicable to all liquid crystals used in biological NMR to date.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
Measurement of global orientational constraints, such as residual dipolar couplings (RDCs), offers many new opportunities in NMR studies of macromolecules in solution, ranging from structure validation and rapid structure determination to the study of dynamics (Clore et al., 1998; Cornilescu et al., 1998; Delaglio et al., 2000; Hansen et al., 1998; Hus et al., 2001; Meiler et al., 2001; Tjandra et al., 1997; Tolman, 2002; Tolman et al., 1995). This technology relies on weak alignment of the macromolecule in solution, usually accomplished by means of a dilute liquid-crystalline medium (Barrientos et al., 2000, 2001; Clore et al., 1998; Fleming et al., 2000; Hansen et al., 1998; Prestegard and Kishore, 2001; Prosser et al., 1998; Ruckert and Otting, 2000; Tjandra and Bax, 1997).

Previously, we and others have shown that alignment in nematic phases of neutral (undoped) bicelles is dominated by steric effects, and that the molecular alignment tensor can be predicted accurately on the basis of the molecule's three-dimensional shape (Almond and Axelsen, 2002; Fernandes et al., 2001; Zweckstetter and Bax, 2000). Predicted alignment tensors for a given structure have been used for the differentiation of monomeric and homodimeric states (Zweckstetter and Bax, 2000), conformational analysis of dynamic systems such as oligosaccharides (Azurmendi and Bush, 2002), refinement of nucleic acid structures (Warren and Moore, 2001), determination of the relative orientation of protein domains (Bewley and Clore, 2000), and validation of structures of protein complexes (Bewley, 2001). The ability to predict dipolar couplings for a given protein structure also provides unique opportunities when attempting to classify protein fold families on the basis of unassigned NMR data, potentially increasing data throughput in structural genomics (Valafar and Prestegard, 2003).

Solute alignment can be altered considerably by adding net charge to bicelles (Ramirez and Bax, 1998). As a direct consequence, the prediction of macromolecular alignment based on steric obstruction fails completely for highly charged liquid-crystalline media (Zweckstetter and Bax, 2000) such as the popular nematic phases of phage particles fd and Pf1 (Clore et al., 1998; Hansen et al., 1998).

Analysis of solute alignment observed in apolar organic liquid-crystalline media shows that steric interactions frequently dominate, but electrostatic forces can be important (Dingemans et al., 2003; Terzis et al., 1996). There is still much debate as to which electrostatic interactions are most important. Some studies indicate that in apolar media it is primarily the electric quadrupole moment of the solute that dominates its alignment (Syvitski and Burnell, 1997, 2000). Other studies argue that the dipole moment plays a major role (Photinos et al., 1992; Photinos and Samulski, 1993). In aqueous solution, full calculation of electrostatic forces between charged macromolecules is a very complex and computationally demanding problem (Davis and McCammon, 1990; Emsley et al., 1991; Honig et al., 1993; Vroege and Lekkerkerker, 1992; Warshel and Aqvist, 1991).

Here, we demonstrate that for aqueous lyotropic liquid-crystalline media, such as that of filamentous phage, a highly oversimplified model, which approximates the electrostatic interaction between a solute and an ordered phage particle as that between the solute surface charges and the electric field of the phage, predicts both the magnitude and orientation of the solute's alignment tensor with reasonable accuracy. Our results demonstrate that in charged liquid-crystalline media electrostatic and steric interactions are the two major contributors to solute ordering. In contrast to the case of steric ordering, we find that depending on the system studied the predicted alignment can be quite sensitive to structural details of the model, in particular the precise location of charged groups on the surface of a protein.

A very recent and independent study by Ferrarini proposes a similar method to predict alignment in charged nematic media and compares experimental dipolar couplings of the first Igg-binding domain of protein G (GB1) with couplings predicted on the basis of its third Igg-binding domain (GB3) (Ferrarini, 2003). The program we describe here and make available predicts rather different alignment for the GB1 and GB3, resulting from their different surface charge distributions, but shows reasonable agreement with experimental data over a wide range of ionic strength. The program also permits evaluation of the contributions from net charge, dipole, and quadrupole moments to ordering, and application to four different proteins and a DNA dodecamer indicates that the relative importance of these terms varies strongly with the system studied.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
Sample preparation
Dipolar coupling data were collected for four small proteins dissolved in dilute suspensions of Pf1 bacteriophage: ubiquitin, the RecA-inactivating protein DinI, and the first and third Igg-binding domains of protein G, GB1 and GB3. Pf1 phage was purchased from ASLA (Riga, Latvia; http://www.asla-biotech.com/asla-phage.htm). ¹³C/¹⁵N/²H-enriched ubiquitin was prepared as described by Sass et al. (1999) and kindly provided to us by S. Grzesiek. ¹⁵N-enriched RecA binding protein DinI and ¹⁵N-enriched GB1 and GB3 were prepared as described previously (Gronenborn et al., 1991; Ramirez et al., 2000). NMR samples contained 1 mM ubiquitin, 0.3 mM DinI, 0.2 mM GB1, and 2 mM GB3 in 90% H₂O/10% D₂O, pH 7. Titration studies were also carried out for a selectively ¹³C/¹⁵N-labeled dodecamer, d(CGCGAATTCGCG)₂ (Tjandra et al., 2000) with the uniformly ¹³C-enriched nucleotides indicated by underlines.

NMR spectroscopy and processing
All NMR experiments were carried out at 25°C, in thin-wall Shigemi microcells (Shigemi, Allison Park, PA) on Bruker DRX600 (ubiquitin, DinI, GB1) and DRX800 (GB3, Dickerson dodecamer) spectrometers (Bruker, Billerica, MA) equipped with three-axis pulsed field-gradient ¹H/¹⁵N/¹³C probeheads. Residual ¹⁵N-¹H dipolar couplings in the four proteins were derived from in-phase/anti-phase (IPAP) [¹⁵N, ¹H]-heteronuclear single quantum coherence (HSQC) experiments (Ottiger et al., 1998). One-bond sugar ¹D_CH RDCs in the Dickerson dodecamer were measured from F₁-coupled ¹H-¹³C HSQC spectra. Spectra were processed and analyzed using NMRPipe/NMRDraw (Delaglio et al., 1995).

Determination of Pf1 concentration
Although reasonable quantitative estimates of the Pf1 concentration can be obtained from its ultraviolet absorbance at 270 nm, using an extinction coefficient of 2.25 cm^–1 mg^–1 ml, a more convenient way to characterize Pf1 concentration simply monitors the lock solvent quadrupole splitting, Q_cc(²H). As incomplete Pf1 ordering in the paranematic phase scales Q_cc(²H) by the same factor as the observed dipolar splittings for the solute (Zweckstetter and Bax, 2001), this makes it straightforward to reference all results reported in this study to a single Pf1 concentration, which was chosen to be 17 mg/ml, corresponding to a Q_cc(²H) splitting of 15 Hz.

Salt titrations
The change in alignment strength and orientation was investigated for a variety of ionic strengths. Titrations were started at very low ionic strength and the ionic strength was continuously increased by addition of 1 M (low ionic strength) or 4 M (high ionic strength) NaCl solutions. In total, ¹⁵N-¹H dipolar couplings were measured in DinI at NaCl concentrations of 30, 50, 110, 210, 315, 400, and 500 mM; in ubiquitin at 20, 50, 105, 240, 340, and 450 mM NaCl; in GB3 at 10, 20, 54, 105, 210, 360, and 500 mM NaCl; and in GB1 at 10, 20, 53, 105, 206, 300, 400, and 500 mM NaCl. To overcome problems associated with very strong ordering at low ionic strength, dipolar coupling measurements in ubiquitin, DinI, and GB3 needed to be performed at dilute Pf1 concentrations.

For DinI and GB3, starting concentrations of 12 mg/ml (quadrupole splitting of the solvent ²H resonance Q_cc(²H) = 10.7 Hz) and 9.3 mg/ml Pf1 (Q_cc(²H) = 8.2 Hz), respectively, were used (Zweckstetter and Bax, 2001). Ubiquitin aligns very strongly and the starting Pf1 concentration had to be reduced to 4 mg/ml, with the result that the titration is already started in the paranematic region. To evaluate ordering in the fully liquid-crystalline state, ubiquitin dipolar couplings were also measured using 20.5 mg/ml Pf1 and 240 mM NaCl.

For titration of the Dickerson dodecamer, the sample was initially washed multiple times with an 8 mM NaCl buffer (pH 6.8). RDCs were measured at five different salt concentrations: 8, 50, 100, 210, and 300 mM NaCl. The starting concentration of Pf1 was 26.2 mg/ml (Q_cc(²H) = 23.2 Hz). During the titration series, this quadrupole splitting decreased to 20 Hz, but again, all data have been rescaled to 17 mg/ml Pf1 concentration. Errors in ¹D_CH measurements are below 1 Hz and therefore have negligible influence on the error in alignment tensor, which is dominated by errors in the atomic coordinates.

Analysis of experimental residual dipolar couplings
Alignment tensors in liquid-crystalline or paranematic Pf1 were obtained by fitting one-bond RDCs to high-resolution NMR (for DinI, ubiquitin, and DNA: mean of 1GHH, Ramirez et al., 2000; 1D3Z, Cornilescu et al., 1998; and 1DUF, Tjandra et al., 2000, respectively) or crystal structures (for GB3: residues 6–61 of 1IGD, Derrick and Wigley, 1994; for GB1: 1PGA, Gallagher et al., 1994), using singular value decomposition (SVD) (Losonczi et al., 1998) as implemented in the program PALES (M. Zweckstetter, unpublished). For x-ray structures hydrogens were added using the program MOLMOL (Koradi et al., 1996). Errors in these alignment tensors are mainly influenced by uncertainties in the reference structures. These errors were estimated using a Monte Carlo structural noise method recently introduced (Zweckstetter and Bax, 2002). DinI aligns very strongly at 30 and 50 mM NaCl, 12 mg/ml Pf1, and only nine ¹⁵N-¹H RDCs could be measured. To avoid large errors in ordering strength introduced by structural noise during SVD, the alignment magnitude at these salt concentrations was determined from its value at 100 mM using the scaling of RDCs when going from 100 to 30 or 50 mM NaCl (neglecting any change in alignment tensor orientation).

Implementation of charge/shape prediction
The charge/shape prediction has been implemented into the dipolar coupling analysis software PALES using the C programming language. To make argument lists more concise, default parameters are used in the absence of additional command-line arguments. Thus, the minimal input to the program consists of a file name containing the structure coordinates, the pH of interest, the concentration of the liquid-crystalline medium, and the ionic strength. In addition, electrostatic potentials can be supplied, user defined charges can be read, the order parameter of the liquid crystal and its average surface charge can be adjusted, and experimental dipolar couplings can be compared to those simulated. A charge/shape simulation using the full-charge model for a protein of 80 residues can be performed in about 1 min on a single processor SGI origin (270 MHz) workstation (Silicon Graphics, Houston, TX). Use of a multipole expansion speeds up the calculation, particularly for larger proteins as only a total of ten charges is taken into account.

Calculation of electrostatic interactions
In our simplified model, proteins and nucleic acids, as charged polymers, are represented by the charges of their ionizable residues. The current studies were performed at pH 7. Therefore, all ionizable residues, except for histidine, are assumed to be fully charged (charge of 1.0 e for Arg and Lys, charge of –1.0 e for Glu and Asp). The N- and C-termini also carry a charge of 1.0 e and –1.0 e, respectively. The protonation state of His is calculated by the program using the Henderson-Hasselbach equation, assuming a reference pKa of 6.3 (charge of 0.33 e at pH 7; Antosiewicz et al., 1994). Charges are distributed evenly over the heavy atoms involved (e.g., both N atoms for Arg, but only N for Lys). Considering more detailed charge models for histidines or other side-chain charges (such as deviations from default pK values) did not consistently improve the quality of our charge/shape prediction, but introduces additional simulation parameters and was therefore not used. Similarly, use of partial charges for the remaining atoms, as defined in CHARMM TOPH19 parameter set, is found to have only a small effect (compared to other factors such as the quality of the reference structure) on predicted alignment tensors. Placing effective charges only at the positions of ionizable residues reduces the total number of charges needed to calculate the electrostatic interactions, and therefore speeds up the computation significantly. For example, for DinI, only 46 instead of 502 charges needed to be used.

Obstruction effects were evaluated using only surface accessible atoms. Surface-accessible atoms were determined using a double cubic lattice method (Eisenhaber and Argos, 1993). For calculation of surface accessibility atom radii were taken from Eisenberg and McLachlan (1986). The solvent radius was set equal to 1.4 Å. To evaluate steric clash between the solute and the liquid crystal particle, radii of solute atoms were set to zero. This allows the solute to access all values of the electrostatic potentials, even the very high values close to the surface of the liquid crystal particle.

For bacteriophage Pf1 a cylinder radius of 3.35 nm, a density of 1.46 g/cm³, a surface charge density of –0.475 e/nm², and a liquid crystal order parameter of 0.9 was used (Zimmermann et al., 1986). The static dielectric constant of water was set to 78.29. Electrostatic potentials for Pf1 bacteriophage were once calculated outside of PALES, and stored as values on a grid with a distance spacing of 0.1 Å between grid points. These potentials are read into PALES and a linear weighting scheme is used to interpolate between grid points. In case of planar liquid crystal particles electrostatic potentials are calculated within PALES from the analytical solution of the nonlinear Poisson-Boltzmann (PB) equation. To reduce calculation times, especially in case of very low liquid crystal concentrations, electrostatic interactions were evaluated until the potential had dropped to 10^–5 (k_BT)/e. At and below potentials of 10^–5 (k_BT)/e the Boltzmann factor p_B is negligible and introduction of this cut-off has no influence on simulated alignment tensors.

We tested the convergence of the results with respect to the spacing of the one-dimensional grid along which the geometric center of the molecule is moved. Increasing the spacing between grid points from 0.2 to 0.5 and to 2.0 Å changes the orientation of the alignment tensor predicted for DinI at 0.20 M NaCl by <3°. Similarly, sampling 362 instead of 122 orientations on the unit sphere and simultaneously increasing the sampling in the third dimension from 18 to 36, i.e., sampling a total of 13,032 instead of 2196 orientations, changes the orientation of the alignment tensor predicted for DinI at 0.20 M NaCl by <2°. Therefore, for all results presented in this study, only 2196 orientations were used.

Analysis of charge/shape-predicted alignment tensors
The accuracy of the magnitude of predicted alignment tensors is assessed using a generalized alignment tensor magnitude,

(1)

where with µ₀ the permittivity of vacuum, h Planck's constant, the magnetogyric ratio, and r_PQ the PQ internuclear distance. In addition, D_a values are used, as visualization of the meaning of ordering parameters is frequently carried out using diagonalized alignment tensors. To evaluate to which degree the orientation of experimental alignment tensors can be predicted, we compare predicted RDCs with experimental ones using Pearson's linear correlation coefficient R_P. This directly addresses the question of how well experimental dipolar couplings can be predicted using our simple model, and avoids errors in experimental alignment tensor orientations introduced by SVD. Finally, the change in experimental alignment tensor orientations with increasing ionic strength is monitored by comparison of RDCs observed at low salt with those measured at later stages of the titration.

Software availability
The program used in the present study, which includes all features of the earlier steric prediction, too, is available upon request from M.Z. (mzwecks{at}gwdg.de) for SGI and Linux machines.

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
Computational approach
Our new method is an extension of the steric obstruction algorithm, described previously (Zweckstetter and Bax, 2000). In the original method, the nematogen is approximated by an infinite wall (bicelles) or infinite cylinder (Pf1 bacteriophage), oriented parallel to the magnetic field (z axis). The center of gravity of the solute is moved on a one-dimensional grid, with a spacing between grid points of 0.2 Å, away from the surface of the liquid crystal model. At each step, a set of 2196 different molecular orientations is sampled. These 2196 orientations are obtained in a two-step procedure. First, the z axis of the molecule samples 122 points on a unit sphere that were determined by a double cubic lattice method (Eisenhaber et al., 1995). This provides a highly uniform sampling of the sphere. The remaining deviation from a completely uniform sampling corresponds to a residual alignment that corresponds to a generalized magnitude, G_Mag, of the alignment tensor (Eq. 1) (Sass et al., 1999) that is smaller than 10^–7. With typical experimental alignment strength of G_Mag 10^–3 (see Table 1) this introduces negligible errors into the calculations. In a second step, the molecule is rotated around the z axis in steps of 20°.

(2)

where _i indicates the angle between the i^th molecular axis and the z axis (magnetic field direction) and _ij the Kronecker delta. In the original, steric model the overall molecular alignment tensor, A^mol, is simply the linear average over all nonexcluded A matrices (Eq. 2). In the electrostatic extension of this model, each nonexcluded A matrix is weighted according to its Boltzmann probability, after calculating the corresponding electrostatic potential of the solute (see below). Using periodic boundary conditions, sampling is restricted to distances r between the solute center of gravity and the center of the bilayer or cylinder for which r < d/(2V_f) (wall model), or r < d/(4V_f)^1/2 (cylinder), where d is either the wall thickness (40 Å for bicelles) or the cylinder diameter (67 Å for Pf1), and V_f is the nematogen volume fraction. The imperfect alignment of liquid crystals is taken into account by multiplication of A^mol with the order parameter of the liquid crystal.

Continuum electrostatic theory (Debye and Hueckel, 1923) is used for calculating the electrostatic interaction energy: the solute is embedded in a dielectric medium containing excess ions, in addition to the counterions neutralizing the solute and nematogen. The nonlinear PB equation is used to derive the electrostatic potential (Chapman, 1913; Gouy, 1910). Even within the simplifications of a continuum description, calculation of the electrostatic potential would require solving a full three-dimensional electrostatics problem for each distance and orientation of the solute with respect to the surface of the charged liquid crystal particle. Instead, we further simplify the problem by treating the solute as a particle in the external field of the liquid crystal. Moreover, we assume the nematogen to carry a uniform charge density instead of discrete surface charges. The nonlinear three-dimensional PB equation is then solved only once, in the absence of the solute, yielding an electrostatic potential (r). The distance and orientation-dependent electrostatic free energy of the protein comprising partial charges q_i at positions r_i is then approximated by

(3)

The Boltzmann factor p_B = exp[–G_el(r,)/k_BT] provides relative electrostatic weights when averaging the individual alignment tensors, derived for each orientation and distance (Eq. 3), to yield an overall solute alignment tensor

(4)

For a flat surface an analytical solution of the nonlinear PB equation exists (Chapman, 1913; Gouy, 1910). For uniformly charged cylinders (bacteriophage) the method of Stigter (1982a,b) is used assuming symmetric monovalent ions and vanishing potential at infinity. We also solved the PB equation with Neumann boundary conditions, i.e., zero normal derivative on an outer cylinder of radius r, to mimic the situation in a two-dimensional, hexagonally packed lattice of phages with neutral unit cells (Hummer et al., 1998). This, however, does not lead to significant changes. For an outer radius of r = 200 Å and salt concentrations above 0.01 M, the differences between the two solutions are confined to the outer boundary region, where the potentials are already very weak (data not shown). This can be understood using charge neutralization arguments: If the Debye length 1/ is much smaller than r, the phage is charge neutralized well before reaching the outer boundary. The two boundary conditions lead to substantially different results only if 1/ is similar to r or larger.

Prediction of charge-induced molecular alignment
Experimental RDCs have been measured in four different proteins and a DNA oligomer d(CGCGAATTCGCG)₂, often referred to as the Dickerson dodecamer (Dickerson and Drew, 1981). For each molecule a high-resolution x-ray or NMR structure is available and experimental RDCs fit well to the structures when evaluated by SVD (Table 1; Fig. 1). These RDCs and structures were used to evaluate the validity of Eq. 4.

View larger version (18K): [in this window] [in a new window]   FIGURE 1  Experimental 1DNH values at 0.20 M NaCl (solid line) and 1DNH values best-fitted to the 1.1 Å crystal structure of GB3 (PDB code: 1IGD) (dashed line). The dipolar coupling quality factor Q (Ottiger and Bax, 1999) is 10%.

View larger version (17K): [in this window] [in a new window]   FIGURE 2  Comparison of experimental 1DNH values with values predicted on the basis of the molecules' three-dimensional shapes (not taking into account electrostatic effects) for ubiquitin (A), DinI (B), GB3 (C), GB1 (D), and the Dickerson dodecamer (E). Experimental values were obtained at 0.20 M NaCl, pH 7, 25°C, using Pf1 bacteriophage as alignment medium. Experimental and predicted 1DNH values were rescaled to 17 mg/ml Pf1. High-resolution NMR (PDB codes for DinI, ubiquitin, and DNA: 1GHH, 1D3Z, and 1DUF, respectively) or crystal structures (for GB3: residues 6–61 of 1IGD; for GB1: 1PGA) were used for shape prediction. Straight lines correspond to y = x.

Does inclusion of electrostatic interactions, treated in the highly simplified manner described above, improve the prediction of residual dipolar couplings from a known three-dimensional structure? The answer to this question is obtained by inspection of Fig. 3. For all four proteins much improved agreement between predicted and observed RDCs is observed compared with the shape-predicted values of Fig. 2. For the 81-residue recA-binding protein DinI, experimental ¹D_NH values agree almost perfectly with those predicted on the basis of DinI's charge/shape model (Ramirez et al., 2000), with a Pearson's correlation coefficient R_P = 0.96 (Fig. 3 B). In addition, the alignment magnitude of G_Mag = 3.2 x 10^–3 predicted at 17 mg/ml Pf1 is close to the experimental value of G_Mag = 4.1 x 10^–3 (Table 1). With R_P = 0.95, the prediction quality for ubiquitin is comparable to that of DinI. The magnitude of ubiquitin's alignment is also much better predicted by the charge/shape model than by the pure shape model (compare Fig. 3 A with Fig. 2 A). For GB3 and GB1 the quality of the predictions is lower, with R_P = 0.88 and R_P = 0.75, respectively.

View larger version (17K): [in this window] [in a new window]   FIGURE 3  Comparison of experimental 1DNH values with values predicted on the basis of the molecules' three-dimensional shapes and charge distributions (taking into account electrostatic effects) for ubiquitin (A), DinI (B), GB3 (C), GB1 (D), and the Dickerson dodecamer (E). For other details, see legend to Fig. 2.

View larger version (12K): [in this window] [in a new window]   FIGURE 4  Correlation between experimental 1DNH values and charge/shape-predicted couplings. Results for five molecules with very different electrostatic properties are shown: (A), ubiquitin; (B), DinI; (C), GB3; (D), GB1; and (E), Dickerson dodecamer. In (F) the charge/shape prediction results for the ensemble of 3GB1 structures is shown. Dashed line indicates RP values for the regularized average structure of the ensemble. Alignment magnitudes are presented in Fig. 7.

View larger version (22K): [in this window] [in a new window]   FIGURE 7  Dependence of alignment magnitude on ionic strength for ubiquitin (A), DinI (B), GB3 (C), GB1 (D), and the Dickerson dodecamer (E). Filled symbols mark values obtained from experimental dipolar couplings by SVD. Error bars represent the influence of structural noise on SVD calculations (Zweckstetter and Bax, 2002). Solid gray lines represent the alignment magnitudes for the full electrostatic model without multipole scaling using the following structures: (A), 1D3Z; (B), 1GHH; (C), residue 6-61 of 1IGD; (D), 1PGA; and (E), 1DUF. Solid black lines correspond to the charge/shape predictions that were obtained with empirically optimized scaling factors for the monopole, dipole, and quadrupole moment (see text). Error estimates (dotted lines) are standard deviations obtained from the ensemble of NMR structures (A,B) and amount to roughly 30%. The alignment magnitude, GMag, is defined in Eq. 1. To facilitate comparison with the familiar Da(NH) values, GMag values have been multiplied by

A first indication of the sensitivity of charge/steric predictions to the conformation of charged side chains is seen when evaluating prediction results for the entire ensemble of 20 NMR structures of DinI deposited in the Protein Data Bank (PDB). Despite the low root mean-square (RMS) spread of the highly refined structure (0.17 Å for the backbone, and 0.84 Å for all atoms) and the high quality of backbone cross-validation and Ramachandran map statistics, a considerable degree of variation in the alignment orientation and magnitude is observed among these 20 members of the ensemble (Fig. 5). In particular, two of the 20 NMR models result in a significantly worse prediction, with R_P = 0.88 (compared to an average value of 0.96). Despite the high quality and small RMS spread of the backbone structure, there is considerable variation in the predicted alignment magnitude, too. It covers a D_a(NH) range from 13 to 26 Hz, with an average value of 19.7 Hz. The lowest RDC correlation factor, R_P = 0.88, is obtained for model 15 of the NMR ensemble (PDB code: 1GHH; Ramirez et al., 2000). It can be traced back to a change in the ₁ angle for Lys⁵⁷. Adjusting the single ₁ angle of Lys⁵⁷ to –60° (as observed in most of the other NMR structures) increases R_P to 0.95.

View larger version (20K): [in this window] [in a new window]   FIGURE 5  Distribution in magnitude, Da(NH), rhombicity, Dr/Da, dipolar coupling correlation, RP, and orientation of charge/shape-predicted alignment tensors, resulting from small structural differences across the NMR ensemble of 20 DinI structures (PDB code: 1GHH) at 0.20 M NaCl, pH 7. The prediction quality is evaluated using Pearson's linear correlation factor, RP. The deviation in orientation is expressed as a generalized angle, according to Sass et al. (1999).

Multipole expansion of protein charge distribution
Results obtained for the highly negatively charged Dickerson dodecamer (net charge –22 e), however, indicate some problems with our simple model. Although the orientation of the alignment tensor is accurately predicted (RDC correlation R_P = 0.95), the alignment magnitude is underestimated by about a factor of 2 (Fig. 3 E). Similarly, dipolar couplings predicted for ubiquitin—although much improved compared to a pure obstruction model—still underestimate experimentally observed values (Fig. 3 A).

To investigate the differences between predicted and experimental alignment magnitudes and to find out which electrostatic interactions are most important for molecular ordering, we dissected the electrostatic contribution to molecular alignment into its multipole components, up to third order. The quadrupole tensor Q = (Q_ß) with , ß = 1, 2, 3 of a set of point charges q_i (i = 1,...,N) at positions r_i = (x_i,y_i,z_i)^T is defined as

(5)

where r_i1 = x_i, r_i2 = y_i, r_i3 = z_i, and _ß; is the Kronecker delta which is 1 for = ß and zero otherwise. Note that the trace of Q is zero, The positions of the charges r_i in Eq. 5 are given with respect to a reference center. The center defines the origin of the Cartesian coordinate system. If the charge distribution has a nonzero net charge or a nonzero dipole moment, then the quadrupole tensor Q is not uniquely defined, i.e., it depends on the choice of the reference center. In all calculations below, the reference center is chosen to coincide with the center of mass.

To calculate the electrostatic interaction energy of a quadrupolar system in an external electrostatic potential (r), we seek to replace the quadrupole moment of all charges q_i with a simpler set of charges that has the same quadrupolar tensor and no net charge and dipole moment. First, we calculate the eigenvalues _x, _y, and _z of Q with corresponding orthonormal eigenvectors v_x, v_y, and v_z, such that

(6)

with

and

(7)

i.e., where = 1, 2, and 3 correspond to x, y, and z. One can then place charges q_x at ± v_x, q_y at ± v_y, q_z at ± v_z and –2(q_x + q_y + q_z) at the origin (i.e., the reference center, for example the center of mass). This neutral arrangement of seven charges has the same quadrupole tensor Q, and a vanishing dipole moment, if

(8)

The electrostatic interaction energy U of the seven charges with the external field (r) is then:

(9)

In the limit of 0, Eq. 9 reduces to an expression in terms of the electric field gradient <EFG>. Here, an appropriate choice for the distance from the center is = 1 Å.

Fig. 6 shows the correlation between experimental RDCs and those predicted based on a model that takes into account all protein/DNA atoms for evaluation of steric obstruction with the liquid crystal particle, but uses a simplified model for the charge distribution of the biomolecule. In this simplified electrostatic model the charge distribution is reduced to only the monopole q (dashed lines), only the dipole µ (dotted lines), only the quadrupole Q (dash-dot-dash lines), or to a combination of all three multipole moments (solid lines). As expected, it is not sufficient to consider only the monopole component. Electrostatic interactions with the directionless monopole affect the orientational distribution of the solute only by reweighting the distance distribution with respect to the charged surface, and thus predominantly influence the alignment magnitude and not the orientation of a molecular alignment tensor. The residual correlation observed in Fig. 6, especially at high salt, is due to the nonzero steric contribution to alignment.

View larger version (17K): [in this window] [in a new window]   FIGURE 6  Correlation between experimental 1DNH values and charge/shape-predicted couplings as a function of ionic strength for increasing complexity of the charge model used for biomolecules. In the simplified electrostatic model the charge distribution is reduced to only the monopole q (dashed lines), only the dipole µ (dotted lines), only the quadrupole Q (dash-dot-dash lines), or to a combination of all three multipole moments (solid lines). Results are shown for four proteins: (A), ubiquitin; (B), DinI; (C), GB3; (D) GB1.

Similar conclusions are obtained from the relative magnitude of alignment that each multipole moment contributes to the total alignment strength. Alignment of biomolecules in charged liquid crystals is generally a function of all three multipole moments. Moreover, their relative importance seems to depend on ionic strength with the quadrupole moment Q becoming more important at higher salt concentrations (data not shown).

Empirical optimization of multipole contributions
The alignment tensor direction remains close to experimentally observed values for all five molecules across the complete salt range from 0.02 to 0.50 M NaCl (Fig. 4). However, as mentioned before, the alignment magnitude of the Dickerson dodecamer is underestimated by about a factor of 2 (Fig. 3 E) using the full charge distribution. In addition, outside of the 0.10 to 0.20 M salt range, alignment magnitudes predicted for DinI, ubiquitin, GB3, and GB1 tend to scale incorrectly with ionic strength, but to different degrees for different proteins. For ubiquitin the full charge model overestimates alignment magnitude at low ionic strength, whereas for the other proteins it is underestimated (solid gray line in Fig. 7). Ubiquitin is the only one of the four proteins that possesses a positive net charge at pH 7 (q = 0.33 e), suggesting that electrostatic attraction and repulsion are overestimated in our calculations. For DinI, on the other hand, the alignment magnitude is underestimated at high salt concentrations, pointing to incorrect weighting of the dipole and quadrupole moments. To evaluate whether these deviations can be corrected for empirically, we introduce scale factors kM, kD, and kQ by which the pseudocharges of the monopole, the dipole, and the quadrupole moment, respectively, are multiplied. Thus, if we want to decrease the influence of the monopole, we have to multiply the pseudocharge that is placed at the geometric center of the molecule (with a value of 0.33 e for ubiquitin) by a factor smaller than 1. As we have three different variables (kM, kD, and kQ) and five different molecules, optimum scale factors were determined with a three-dimensional grid search. kM, kD, and kQ were varied independently from 1/k to k with k = 3 and a step size for k of 0.2. For each of the 9261 {kM, kD, kQ} combinations, charge/shape predictions are made for all five molecules. All {kM, kD, kQ}; combinations were ranked according to the overall, normalized RMSD between experimental and charge/shape-predicted RDCs, calculated for the five molecules. The {kM, kD, kQ{ combination with the lowest overall RMSD (sum of individual protein RMSDs normalized by each protein's magnitude of alignment) is therefore the one that best reproduces both the magnitude and the orientation of the alignment tensor for all five molecules. According to Fig. 6 the scale factors will also depend on the ionic strength. Therefore, the grid search is repeated for seven salt concentrations, ranging from 0.02 to 0.50 M NaCl. The weight factors that optimize agreement between observed and predicted dipolar couplings are listed in Table 3.

View larger version (16K): [in this window] [in a new window]   FIGURE 8  Correlation between experimental 1DNH values and charge/shape-predicted couplings (solid lines), or SVD best-fit couplings (dotted lines) as a function of ionic strength. The dashed lines indicate the correlations between experimental couplings at a certain ionic strength and those observed at low salt (0.02 M NaCl). Results for five molecules with very different electrostatic properties are shown: (A), ubiquitin; (B), DinI; (C), GB3; (D), GB1; and (E), Dickerson dodecamer. Note that for ubiquitin, at the low Pf1 concentrations used, the drop in RP above 0.30 M NaCl is dominated by the reduced relative accuracy of dipolar couplings. Note that these results were obtained with empirically optimized scaling factors for the monopole, dipole, and quadrupole moment (see text).

View larger version (14K): [in this window] [in a new window]   FIGURE 9  (A) Magnitude of DNA dodecamer alignment relative to the total alignment magnitude, as a function of distance between the center of gravity of the dodecamer and the axis of the viral rod, calculated using shells of 5 Å thickness for various ionic strengths. The x axis indicates the distance from the center of a Pf1 rod to the inner surface of the 5 Å shell. Different NaCl concentrations are indicated in the following way: solid line, 0.01 M; dashed line, 0.1 M; and dotted line, 0.3 M. In gray, the shell distribution for the steric case is shown. (B) Electrostatic potentials of bacteriophage Pf1, using a cylinder radius of 3.35 nm and a surface charge density of –0.475 e/nm2, calculated from the nonlinear Poisson-Boltzmann equation at 0.01 M (solid line), 0.1 M (dashed line), and 0.3 M (dotted line) NaCl.

GB1 (q = –4.0 e; µ 80 Debye) shows a much weaker dependence of alignment strength on salt concentration than DinI and GB3. Although it differs from GB3 by only three charged residues the alignment magnitude decreases by just 40% when going from 0.02 to 0.50 M NaCl (Fig. 7 D). In addition, for salt concentrations above 0.15 M NaCl the steric component of the alignment tensor plays a significant role. This causes an experimentally observed rotation of the alignment tensor orientation by 60° when going from 0.02 to 0.50 M NaCl. Accordingly, experimental dipolar couplings measured at 0.02 M NaCl and at 0.50 M NaCl correlate poorly with one another (R_P = 0.44). The quality of charge/shape-predicted dipolar couplings in GB1, on the other hand, remains high across the entire salt range and even increases at high ionic strength (Fig. 8 D). This demonstrates that the increasing steric contribution to alignment is correctly taken into account by our simple model. As discussed earlier, the overall lower correlation for the protein G domains compared to ubiquitin and DinI can be attributed largely to differences in side-chain orientation when using a static x-ray model to describe the solution structure.

Figs. 7 E and 8 E show the results of charge/shape prediction for a very highly charged molecule, the Dickerson dodecamer (q = –22 e). As discussed earlier, the orientation of the simulated alignment tensor remains unchanged when going from 0.01 to 0.30 M NaCl, in agreement with experimental results. In addition, experimental and simulated RDCs are highly correlated, with a Pearson's correlation factor of 0.96, demonstrating that the orientation of the alignment tensor is correctly predicted.

Remarkably, both our model and the experimental results indicate very little dependence of alignment magnitude on ionic strength. In contrast to the four proteins, the alignment magnitude of the DNA does not decrease with increasing salt concentration. It actually shows a small increase (15%) when going from 0.01 to 0.30 M NaCl (Fig. 7 E). This increase in alignment magnitude can be rationalized by evaluating how much of the total alignment is contributed by molecules as a function of distance from the center of the Pf1 rod. To this extent, we calculated the percentage of total alignment that comes from molecules for which the center of gravity is located within shells of 5 Å thickness, for example within distances ranging from 60 to 65 Å from the Pf1 rod (Fig. 9). No alignment comes from positions very close to the Pf1 particle, as all possible orientations are sterically forbidden. As soon as the molecule is far enough away and orientations exist that are not sterically forbidden, the contribution to overall alignment quickly rises, reaches a maximum, and then drops to zero. If only steric effects are considered, the contributing region is restricted to translational positions where obstruction between the molecule and the hard-wall potential can occur. Outside this region all orientations are equally probable and no preferred orientation exists (shown in gray). When electrostatic interactions are taken into account, the free energy of the protein is orientation-dependent and the main contribution to alignment comes from outside the steric overlap region. The position of the maximum depends on ionic strength: the lower the salt concentration the further away is the shell from the Pf1 rod that contributes the most to ordering. Distances closer to the phage are unfavorable as electrostatic repulsion between the highly negatively charged DNA and Pf1 dominates. Thus, the increase in alignment before it drops to zero results from the initial increase in the number of molecules per shell as a function of the distance from the Pf1 surface, whereas the drop-off in G_Mag per shell initially is small. The same argument explains the Dickerson dodecamer's increase in alignment magnitude with increasing salt concentration: at higher ionic strength, more molecules are close to the Pf1 particle, enabling a bigger contribution to alignment by steric obstruction.

	DISCUSSION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
Electrostatic model
Our model only takes into account short-range steric and long-range electrostatic interactions. Nonpolar interactions, such as a Born repulsion and a nonpolar attraction, are neglected. This simplification was motivated by the fact that previous studies have shown that the two opposing nonpolar interactions approximately balance each other (BenTal et al., 1997). In addition, whereas the long-range electrostatic interaction does not significantly depend on the atomic details of the liquid crystal surface, both nonpolar interactions—as they become important when the solute approaches the liquid crystal particle—strongly depend on the model used for the liquid crystal surface. Thus, to include nonpolar interactions, the liquid crystal particle would need to be represented in atomic detail. The good agreement between experimental and predicted alignment tensor values obtained in our study demonstrates that such nonpolar effects at most play a minor role.

There are at least four other simplifications that enter our treatment of electrostatic effects:

1. The molecular solvent—water and ions—is treated at the continuum level, thus excluding the possibility of specific molecular interactions such as ion binding to the protein or phage surface. Although simulation studies of simple models of ions near polyelectrolytes have shown reasonable agreement with such continuum descriptions (Mills et al., 1985; Murthy et al., 1985), we expect deviations in the limits of (a), close contact between protein and interface; (b), low salt concentrations resulting in tight ionic interactions; and (c), high salt concentrations where ion packing contributions play a role. The effects of multivalent and molecular ions are also expected to be described relatively poorly by the simple PB continuum model.
   
2. The orientation-dependent interactions of the protein with the charged interface are calculated by using the electrostatic potential of the unperturbed interface, without calculating the combined effect of protein and interface on the distribution of ions in their vicinity.
   
3. Although our model explicitly treats the shape and charge distribution of the solute, the charge distribution on the liquid crystal particle is assumed to be uniform. This approximation is expected to lead to errors when the protein is close to the interface, where specific interactions with interfacial charge distribution are important.
   
4. The protein structure and the charge distribution on the protein and interface are assumed to be rigid and independent of orientation and distance. Errors are expected for proteins with flexible charged side chains that can adjust their conformation to optimize electrostatic interactions with the interface.
   However, the reasonably quantitative agreement between observed and predicted alignment suggests that the above simplifications do not have dramatic adverse effects on the values predicted by our simplified model.
   

The effect of transient binding of the solute to the liquid crystal particle is negligible. If such transient binding occurred, it could introduce important contributions from higher multipole terms as well as nonlinear interaction of separate steric and electrostatic contributions. However, from an NMR perspective, the liquid crystal particle behaves like a solid, and transient binding to it would dramatically increase its transverse relaxation rates, contrary to what is seen for the majority of proteins dissolved in such systems. Even in cases where apparent line broadening in Pf1 medium is observed (DinI and ubiquitin <0.1 M NaCl), such line broadening frequently results from large residual dipolar couplings, whereas ¹⁵N relaxation rates are nearly unchanged compared to isotropic solutions. However, at low ionic strength an increase in the ¹⁵N transverse relaxation rate is observed for ubiquitin (data not shown), together with very strong alignment of the protein, indicating that in this case transient binding events are no longer negligible.

Surface charge of liquid crystal particles
For applying our charge/shape model to structure determination of biomolecules, a potential problem is that the surface charge for many of the liquid crystals is not known accurately and depends on the method of preparation. For example, the charge on the surface of the cellulose-based phase (Fleming et al., 2000) depends on particle size and on the length of acid hydrolysis used. In doped bicelles, protein alignment data indicate a surface charge that is nearly threefold lower than expected on the basis of the mole fraction of added charged amphiphiles, when using an infinite bilayer as the model for the liquid-crystalline medium. The smaller effective charge may result from a preferential partitioning of the charges near or in the rim of the pores, where their effect on solute alignment is very different from those partitioned in the regular bilayer. The same is found for the ostensibly neutral undoped bicelles, where empirically a very small amount of negative charge is found to improve the agreement between predicted and observed dipolar couplings over the purely steric alignment model. The presence of such charge had been predicted (Losonczi and Prestegard, 1998), but the precise magnitude is found to depend on the source of the materials used in the preparation of the liquid crystal. Knowledge of a rather accurate value of the surface charge of the liquid crystal will be important for applications that employ only slightly charged liquid crystals. For such liquid crystals the exact balance between obstruction and electrostatic interaction will determine the orientation and magnitude of the alignment tensor. For highly charged liquid-crystalline media such as bacteriophage Pf1, purple membrane fragments, and cetylpiridinium-based media, on the other hand, exact knowledge of the surface charge of the liquid crystal poses less of a problem. For salt concentrations up to 0.10–0.20 M NaCl, electrostatic interactions will dominate the molecular ordering for most molecules.

Accuracy of charge/shape-predicted alignment tensors
The full charge model, i.e., the model where the protein/nucleic acid is represented by effective charges placed at the positions of its ionizable residues, is superior to the multipole expansion model for predicting the alignment tensor orientation. This is expected, as a multipole decomposition is only valid for distances between the liquid crystal particle and the molecule that are large compared to the spatial extension of the molecule's charge distribution. The magnitude of alignment, on the other hand, is only correctly predicted after replacing the effective charges by a multipole expansion and introduction of empirical scale factors for the monopole, dipole, and quadrupole moment. These empirical scale factors were determined from experimental data obtained for five biomolecules with very different electrostatic properties, and are therefore expected to be applicable to a broad range of proteins and nucleic acids. Nevertheless, the scaling factors must be seen mainly as a compensation tool for the simple nature of our model. The choice of basing these scaling factors on the mono-, di-, and quadrupole moments was made on the basis of their different potential versus distance relationships when varying ionic strength. It cannot be excluded, however, that other empirical parameters could yield equally satisfactory or even better results. Further improvements over our model would require solving the nonlinear PB equation for every distance and orientation of the protein with respect to the surface of the liquid crystal particle.

More importantly, our results indicate that—from the point of application toward structure determination of biomolecules—the accuracy of prediction is ultimately limited by factors other than the details of the electrostatic model. As discussed above, minor rearrangement of the side chains of charged surface residues can dramatically affect the correlation between observed and predicted dipolar couplings. This effect was particularly dramatic for GB1, which has a relatively uniform charge distribution based on its x-ray structure that results in small dipole and quadrupole moments. On the other hand, variation of simulation parameters such as the resolution of the orientational or translational grid, the boundary condition used for calculating the electrostatic potentials, or the detail of the charge distribution, i.e., putting charges only on the side chains of ionizable residues versus using the CHARMM TOPH19 parameter set which puts partial charges on all atoms, has only a small effect on the results obtained from charge/shape prediction. We thus expect that using a more sophisticated and computationally more costly electrostatic model (such as the one implemented in DELPHI; BenTal et al., 1997, 1996; Gallagher and Sharp, 1998; Murray and Honig, 2002; Peitzsch et al., 1995) may only marginally improve the accuracy of predicted alignment tensors—at least not with the quality of three-dimensional structures and their conformational dynamics currently available. NMR relaxation experiments indicate that many of the long side chains of charged residues such as Lys, Arg, and Glu are frequently subject to conformational averaging. In addition, long, flexible, charged N- or C-termini of protein backbones pose a problem for accurate prediction of charge-induced alignment. In principle, it may be possible to use the agreement between observed and charge/shape-based alignment tensor as a restraint when calculating the ensemble, or as a selection criterion when evaluating which ensemble members are acceptable.

From the presented results it is clear that the quality of prediction varies between proteins. This is due to either the oversimplified nature of our model or inaccuracies in the structural model. Therefore, it is difficult to quantitatively evaluate how good a structure is with the simple model introduced here. However, most applications of the charge/steric prediction model (similar to the already widely used steric prediction model) will focus on a distinction between several structural models that are not easily distinguished by other methods. For example, if very different relative conformations (e.g., parallel versus antiparallel) of a homodimeric system are in question, the models are likely to yield very different R_P values. Note that even in the very worst of our cases, outside the limits where reliable data could be expected, R_P > 0.5, whereas in the absence of structural information, chances are 50% that R_P < 0. For example, if for model (A) R_P = 0.8 and for model (B) R_P = –0.3 then the protein is very likely of the model (A) type.

Interactions responsible for molecular alignment of biomolecules
There has been much debate about which electrostatic interactions are most important for solute alignment in apolar liquid-crystalline media (Dingemans et al., 2003; Emsley et al., 1991; Photinos et al., 1992; Photinos and Samulski 1993; Syvitski and Burnell, 1997, 2000; Terzis et al., 1996; Vroege and Lekkerkerker, 1992). In aqueous media the situation is fundamentally different as the ion distribution in the solvent dramatically perturbs the electrostatic potential.

For electrostatic interactions, higher multipole orders become