Molecular Dynamics Simulations of Aqueous Pullulan Oligomers

Feb 22, 2005 - on aqueous solutions of four oligomeric segments of the glucan pullulan: the trimer G3 (comprising one polymer repeating unit), the hex...
0 downloads 0 Views 616KB Size
Biomacromolecules 2005, 6, 1239-1251

1239

Molecular Dynamics Simulations of Aqueous Pullulan Oligomers Simon Jaud, Douglas J. Tobias, and David A. Brant* Department of Chemistry, University of California, Irvine, California 92697-2025 Received September 1, 2004; Revised Manuscript Received January 18, 2005

Molecular dynamics (MD) simulations have been used to model small-angle X-ray scattering (SAXS) data on aqueous solutions of four oligomeric segments of the glucan pullulan: the trimer G3 (comprising one polymer repeating unit), the hexamer (G3)2, the nonamer (G3)3, and the dodecamer (G3)4. The AMBER* force field was used in conjunction with the GB/SA continuum solvation model to calculate both the mean global dimensions of the oligomers from the limiting small angle scattering behavior and the shorter range structural information implicit in the Debye scattering function at larger scattering angles. This same force field and solvation treatment were employed earlier by Liu et al. (Macromolecules 1999, 32, 8611-8620) with apparent success in a rotational isomeric state (RIS) treatment of the same experimental data. The present work discloses that, despite numerical success in modeling the SAXS data, the RIS treatment, which includes only the interactions within dimeric segments of the polymer chain, fails to account accurately for excluded volume effects at the range of 3-12 sugar residues in the polymer backbone. It is suggested that MD simulations using continuum solvation models can be used to circumvent errors inherent in the computationally efficient RIS treatments of polymer nano- and picosecond dynamics while at the same time avoiding the heavy computational requirements of all-atom methods. Introduction Pullulan is a water soluble, random coil glucan that serves as a paradigm for the behavior of aqueous polysaccharides.1-3 It is a regularly repeating copolymer with the chemical structure {f6)-R-D-glucopyranosyl-(1f4)-R-D-glucopyranosyl-(1f4)-R-D-glucopyranosyl-(1f}n.4 A hexameric segment showing a single central (1f6)-linkage is shown in Figure 1. The polymer may thus be viewed as a succession of R-(1f6)-linked (1f4)-R-D-triglucosides, i.e., maltotriose units (G3). We have previously investigated the properties of aqueous pullulan using viscometry and light scattering5 and small-angle X-ray scattering (SAXS)6 to provide the basis for developing a detailed connection between the chemical structure of pullulan and its physical behavior in aqueous solution. Molecular models have been developed at the atomistic level to interpret the experimental observations.6-10 A recent rotational isomeric state (RIS) treatment accounting only for nearest neighbor interactions along the chain backbone successfully reproduced experimental SAXS data on the four pullulan oligomers G3, (G3)2, (G3)3, and (G3)4, containing respectively 3, 6, 9, and 12 sugar residues; the same model matched the measured unperturbed dimensions of high molecular weight pullulan.6 This model was based on an estimation of the conformational energy surface of representative dimeric segments of the pullulan chain using the AMBER* molecular mechanics force field tuned for application to carbohydrates; aqueous solvation was included * To whom correspondence should be addressed. Telephone: +1 (949) 824-7840. Fax: +1 (949) 824-9181. E-mail: [email protected].

Figure 1. A hexameric segment of pullulan showing the (1f6)- and (1f4)-glycosidic linkages. Only carbon and oxygen atoms are shown.

using the generalized Born/surface area continuum treatment (GB/SA).11,12 This AMBER* + GB/SA method was implemented within the commercial MacroModel 5.5 software package.13 More detailed RIS models of the pullulan oligomers G3 and (G3)2 which took into account longer range interactions within the oligomers were based on the CFF95 force field in conjunction with the commercial CERIUS2 software package.10 This treatment was coupled with a Poisson-Boltzmann (PBF) model for aqueous solvation14,15 available in the commercial JAGUAR molecular modeling package.16 This approach was somewhat less successful than the nearest neighbor model in matching the SAXS data for the hexamer (G3)2.10 In the present work, we use molecular dynamics (MD) simulations to explore the aqueous conformational properties of the oligomers G3, (G3)2, (G3)3, and (G3)4 at 300 K by employing the AMBER* + GB/SA force field implemented within the MacroModel 7.0 package. This approach permits exploration of the role of excluded volume effects on pullulan

10.1021/bm049463d CCC: $30.25 © 2005 American Chemical Society Published on Web 02/22/2005

1240

Biomacromolecules, Vol. 6, No. 3, 2005

Figure 2. Disaccharides R-maltose and R-isomaltose. The glycosidic torsion angles are defined as φ4 (O5-C1-O1-C4′) and ψ4 (C1-O1C4′-C5′) for the (1f4)-linked maltose and φ6 (O5-C1-O1-C6′), ψ6(C1-O1-C6′-C5′), and ω6 (O1-C6′-C5′-O5′) for the (1f6)-linked isomaltose. The three exo angles are τ5 (O5-C5-C6-O6) and τ5′ (O5′-C5′-C6′-O6′) for maltose and τ5 (O5-C5-C6-O6) for isomaltose. The staggered conformations of τ5 and τ5′ are traditionally designated gg (-60°), gt (60°), and tg (180°). The H6 and H6′ atoms have been omitted for clarity.

oligomer conformation using a force field that provided excellent correlation of experimental SAXS and light scattering results on pullulan oligomers and the pullulan high polymer in the RIS approximation.6 The RIS method ignores intrapolymeric interactions beyond the range of first neighboring sugar residues and is thus applicable only to relatively short polymer chains or to measurements made on unperturbed high polymers in a Flory θ solvent. Moreover, MD simulations with AMBER* + GB/SA afford an opportunity to test the applicability of MD with continuum solvation with its inherent gain in computational speed and its ability to explore processes with longer relaxation times than typically accessible with MD.17 For comparison, we also carried out MD simulations using a CHARMM force field designed for carbohydrates18,19 with explicit TIP3P20 aqueous solvation. Details are not reported here. Pullulan Structural Geometry. The characteristic dimeric segments of pullulan are [fx)-R-D-glucopyranosyl-(1f4)R-D-glucopyranosyl-(1f] and [f4)-R-D-glucopyranosyl(1f6)-R-D-glucopyranosyl-(1f], where x may be either 4 or 6 for the (1f4)-linked segment. These segments are shown in Figure 2, where both types of (1f4)-linked segments are represented by the disaccharide R-D-maltose (R-D-glucopyranosyl-(1f4)-R-D-glucopyranose) and the (1f6)-linked segment is represented by the disaccharide R-Disomaltose (R-D-glucopyranosyl-(1f6)-R-D-glucopyranose). Glucose residues are assigned to the preferred 4C1 ring conformation.21 Glycosidic linkage torsion angles φ4, ψ4 and φ6 ,ψ6, ω6 are defined in the legend to Figure 2. The zero

Jaud et al.

and positive direction of rotation conform to standard usage.22-25 Maltose has two exocyclic C5-C6 torsion angles, here designated τ5 and τ5′. Features of the reducing residue of each disaccharide (i.e., the residue containing a free OH1) are designated with a prime. One exocyclic C5-C6 torsion angle is present in isomaltose, designated τ5, and the other C5-C6 torsion angle in isomaltose is the glycosidic linkage dihedral ω6. The staggered conformations of τ5 are traditionally designated gg, gt, and tg. The angle τ5 and its staggered conformers are defined in the legend to Figure 2. We will use the same terminology to describe the staggered conformational states of ω6. There is extensive experimental evidence to show that for R-D-glucans the tg conformer is sparsely populated for both τ526-28 and ω6.29 This is a manifestation of the more general gauche effect.30 A distinct preference of the glycosidic torsion angles φ4 and φ6 of R-D-glucans for values near +60° is a manifestation of the exo-anomeric effect.21,31 Earlier studies have shown that the secondary hydroxyl groups of the glucose residue adopt a “crown” orientation in gas-phase simulations.32,33 This terminology refers to the tendency for all secondary OH bond vectors to be aligned either with the direction of increasing carbon index in the sugar ring (clockwise crown) or against this direction (reverse clockwise crown). Hydrogen bonded (i.e., electrostatic) interactions between adjacent secondary hydroxyls in the gas phase drive formation of the crown arrangements, which in the presence of aqueous solvation may be less strongly preferred.32 The AMBER* simulations with implicit GB/SA water were carried out using the MacroModel 7.0 program. The default cutoff distances were used (van der Waals, 8 Å; electrostatic, 20 Å; hydrogen-bonding, 4 Å), and the dielectric constant was set to 78. The integration time step was set to 1.5 fs, and SHAKE was used to freeze the fast bonds between hydrogens and heavier atoms.34,35 Simulations were done using two 400 MHz processors in parallel on an SGI Octane2/V6. Observable Properties Monitored. The properties monitored during the simulations are those previously measured in SAXS experiments on the oligomers G3-(G3)4.6 In particular, we compute the measurable Debye scattering function P(q) ) I(q)/I(0), where I(q) is the absolute intensity of scattered radiation corresponding to a scattering vector of length q ) (4π/λ) sin(θ/2) and I(0) is the scattering intensity extrapolated to q ) θ ) 0. For any set of atomic scattering centers, I(q) may be calculated from the interatomic distances rij by the Debye formula36 N

I(q) )

∑ ∑ i)1 j)1

〈 〉 sin (rijq)

N

gigjφiφj

rijq

(1)

Here N is the total number of atoms and gi is the net scattering power of atom i, conventionally taken as the atomic number, i.e., the number of electrons of atom i less the mean atomic number of the corresponding volume of background solvent. For water the mean atomic number is 0.336 electrons/Å3 × Vi Å3, where Vi is the volume of atom i computed on the assumption that atom i is a sphere of radius

Biomacromolecules, Vol. 6, No. 3, 2005 1241

Simulation of Pullulan Oligomers

Figure 3. Scattering form factors of carbon (plain line) and oxygen (dashed line) calculated using eqn. 2 with Ri ) fiRiv. The dots correspond to the atomic form factors values given by crystallographic tables.37

Ri. The quantity φi(q) is the atomic form factor for atom i, often approximated by φi(q) )

3[sin(Riq) - (Riq) cos(Riq)] (Riq)3

(2)

This expression for φi(q) assumes uniform electron density within a sphere of radius Ri, which is usually approximated as the van der Waals radius, RiV, of the atom. More accurate values of φi(q) can be obtained from tabulations based on quantum mechanical calculations of atomic electron distribution.37 Following Liu et al.,6 the tabulated form factors for carbon and oxygen were approximated very closely by substituting Ri ) fiRiV in eq 2. For carbon and oxygen, fi ) 0.626 and 0.536, respectively. We used the respective conventional values of RiV for carbon and oxygen, 1.67 and 1.50 Å. The fit to the tabulated atomic scattering factors achieved in this way is shown in Figure 3. Only carbon and oxygen atoms are included in the sum in eq 1; the small contribution of hydrogen atoms to the X-ray scattering is ignored. Note that we also used Ri ) fiRiV to compute Vi in evaluating the mean atomic number of water. The rationale for using fiRiV to compute the volume of water displaced by a given atom is that the atoms displacing water are all bonded, and the nonbonded van der Waals radius overestimates their effective radii. The choice Ri ) fiRiV means that both factors in the products giφi occurring in eq 1 are increased when we take Ri ) fiRiV instead of Ri ) RiV. The experimentally observable mean square radius of gyration can be calculated from the atomic coordinates using eq 3 N

〈s2〉 ) 〈

N

misi2/∑mi〉 ∑ i)1 i)1

(3)

Here si is the distance of atom i from the center of mass of a collection of N atoms and mi is the mass of atom i. For any given conformation of the chain, s2 is the mass weighted, mean square distance of the scattering centers from the center of mass. The angle brackets imply here, and in eq 1, the ensemble average over all orientations and conformations

adopted by the molecules of a macroscopic pullulan oligomer sample. In the context of an MD simulation of a single molecule this ensemble average is equated with the time average of the relevant quantity [s2 or sin(rijq)/rijq] over a time span long enough to sample the configuration space of the molecule ergodically in the equilibrium phase of the simulation. An initial period of relaxation into the equilibrium state (0-10 ns) is excluded from the sampling. Here we will use the notation R ) (s2)1/2 to denote the instantaneous value of the radius of gyration of the subject molecule at any point in the simulation, and Rg(t) will denote the root-mean-square running time average 〈s2〉1/2 after any time t. The long time value of this computed property, Rg ) Rg(t ) ∞), will be compared with the experimental result. An important criterion for determining the time corresponding to t ) ∞ is that Rg(t) has become effectively independent of t. Experimentally one extracts 〈s2〉 from the measured P(q) using the Guinier approximation to P(q)38 P(q) = exp(-〈s2〉q2/3)

(4)

from the limiting slope at q ) 0 of a plot of ln P(q) against q2. The experimental root-mean-square radius of gyration will be designated REg . Because of the inherent approximation in REg , it may be preferable to compare this quantity with the radius of gyration extracted from the simulation by applying the Guinier approximation to the running average of P(q) at t ) ∞ calculated according to eq 1. We designate the theoretical radius of gyration calculated in this way as RGg . The Debye scattering function is frequently analyzed in terms of a plot of q2P(q) vs q, as originally suggested by Kratky.36 This rendition of the data gives increased emphasis to the data at larger q corresponding to interferences arising from shorter range features of the subject molecule. It is often more diagnostic for the specific structural features of the molecule than is P(q) itself. Our experimental data cover the range of q from 0.025 to 0.500 Å-1 and covering the corresponding range of Bragg spacings 250-13 Å.6 Conformational Energy Maps for Dimeric Chain Segments. Relaxed conformational energy maps were prepared for the two characteristic dimeric segments of pullulan, R-Dmaltose and R-D-isomaltose. These maps serve respectively to define the conformational distributions of the (1f4)- and (1f6)-linkages in the absence of interactions beyond the range of nearest neighbor sugar residues and provide a basis for describing the perturbations to these distributions arising from longer range interactions within the oligomers. Results were plotted as energy contour diagrams in φ4, ψ4 space for maltose. For this purpose, the variables φ4, ψ4 were artificially constrained to a regular grid with points located at 10° intervals in the torsion angle space, and the energy was minimized with respect to all other coordinates of the configuration space of the disaccharide. For isomaltose, contour plots were made on a 10° grid in ψ6, ω6, in φ6, ω6, and in φ6, ψ6 space. In each case, the energy was minimized with respect to the other coordinates including the third torsion angle of the (1f6)-linkage, so that the third glycosidic torsional dimension is projected into the plane of the two grid angles. Minimizations with the AMBER* force field

1242

Biomacromolecules, Vol. 6, No. 3, 2005

were done using the Polak-Ribiere conjugate gradient method provided in the MacroModel 7.0 software package. The gradient convergence criterion was set to 0.05 kJ/Å mol, and the grid angles were constrained by a force constant of 100 kJ/Å mol. The disaccharide maps included contributions from the intramolecular terms in the potential energy functions; as specified below, solvation terms in the force field were also included. Starting Conformations and Force Field Parameters. Glucose residues were configured in the preferred 4C1 chair form.21 Atomic coordinates for the glucose residues were drawn from the MacroModel 7.0 structure library, and the residues were subjected to an initial brief relaxation process to remove stress introduced by adoption of these coordinates in conjunction with an arbitrary force field. In all cases, the conformations of the secondary hydroxyl exocyclic Ci-Oi torsion angles were set initially in the crown reverseclockwise orientation suggested by Ha et al.32 to be least energetic for R-D-glucose in the absence of solvent. The C5s C6 torsion angles were set initially to each of the nominal staggered conformations. Nine such combinations τ5-τ5′ exist for maltose, i.e., gg-gg′, gg-gt′, gg-tg′, gt-gg′, gt-gt′, gt-tg′, tg-gg′, tg-gt′, and tg-tg′, whereas for isomaltose there are just three such states of τ5, i.e., gg, gt, and tg. The initial conformations of the C5-C6 torsion angles were retained throughout the minimization process, as was the initial crown orientation. An energy minimized disaccharide map was computed for each of the possible C5-C6 torsional states. A composite disaccharide map was then constructed by assigning to each maltose grid point in the driven torsion angle space, i.e., φ4, ψ4, the energy from that one of the nine minimized maltose maps with the smallest energy at that grid point. Similarly, to each grid point in the driven torsion angle space of the isomaltose (1f6)-linkage was assigned the energy from that one of the three minimized isomaltose maps with the smallest energy at that grid point. It was not deemed necessary for present purposes to spend the effort to obtain truly adiabatic disaccharide energy maps.33 It was evident during preparation of the composite disaccharide maps that for the AMBER* + GB/SA force field the energetic order of the C5-C6 torsional states was gt e gg < tg, consistent with the gauche effect and as expected from earlier experience with this force field.6 These attributes are best displayed by following the τ5 and τ5′ torsion angles during the MD simulations described below. The sugar ring conformations were monitored during the grid dihedral driver process to avoid unwanted transitions in glucose ring conformation away from the preferred 4C1 ring.33 Simulations Using AMBER* + GB/SA Solvation. R-DMaltose and R-D-Isomaltose. Figures 4 and 5 demonstrate that 48 ns MD simulations of the disaccharides R-D-maltose and R-D-isomaltose in water using AMBER* + GB/SA solvation produce trajectories consistent with the composite relaxed disaccharide maps computed with AMBER* + GB/ SA as described above. Minima on the energy surfaces are labeled A4, B4, A6, B6, etc., following the scheme used earlier.6,10 Note that the minima on the isomaltose surface are defined specifically for the ψ6, ω6 projection, which

Jaud et al.

Figure 4. Trajectories for aqueous R-maltose starting from A4 (AMBER*). Contour lines represent increments of 1 kcal/mol.

displays the richest array of minima. Whereas φ6 is the most constrained of the (1f6)-linkage torsion angles and is confined to the approximate range 30 e φ6 e 120° (Figure 5, parts b and c), ψ6 and ω6 explore most of the nine possible staggered conformations for this pair of sequential torsion angles. In the earlier RIS treatments of pullulan,6,10 one rotational isomeric state was associated with each of the minima A6-F6 of the isomaltose surface, whereas a single rotational isomeric state was assigned to minimum A4 of the maltose surface. Hence, in the RIS approximation, all of the tortuosity of the pullulan chain is attributed to the (1f6)linkage, a limitation that is avoided in the present MD approach, even though Figure 4 confirms that neglecting minimum B4 of the (1f4)-linkage is justified with this force field. The maltose surface (Figure 4) is distorted in the vicinity of ψ4 ) 270°, because the dihedral driver routine forced at least one of the glucose residues into a ring conformation other than 4C1 in this domain. Consequently the computed energies for ψ4 > 270° are high and do not accurately describe the behavior of the molecule during the simulation. It is clear that the maltose trajectory penetrates occasionally into this domain of conformation space, and we have verified that this occurs without any change in glucose ring conformation. On the other hand, the reducing residue of isomaltose underwent two ring flips of 1-2 ns duration during the 48 ns simulation. Figure 6 shows the history of the exocyclic torsion angles τ5 and τ5′ for maltose during the trajectory plotted in Figure 4. Clearly the tg state (τ5, τ5′ = 180°) is rarely visited, whereas the gt state (τ5, τ5′ = 60°) is preferred over gg (τ5, τ5′ = -60°) by this force field by a somewhat greater margin than is suggested by experimental observations.26-28 Distributions of the exocyclic C5sC6 torsion angle for the primary hydroxyl for isomaltose and for all of the pullulan oligomers are essentially the same as those shown in Figure 6 when this force field was used. The contour plots in Figures 4 and 5, which are clearly approximate, are

Simulation of Pullulan Oligomers

Biomacromolecules, Vol. 6, No. 3, 2005 1243

Figure 6. Exocyclic torsion angle history for aqueous R-maltose starting from A4 (AMBER*).

Figure 5. Trajectories for aqueous R-isomaltose starting from A6 and τ5 ) gg (AMBER*). Contour lines represent increments of 1 kcal/ mol.

introduced for purposes of illustration only and are not used for computation of observable properties of the oligomers. Pullulan Oligomer G3. Figure 7 reports the results of a 100 ns simulation of aqueous maltotriose G3, the fundamental repeating unit of pullulan. Panel a shows the instantaneous value of R ) (s2)1/2 calculated using eq 3 as a function to the simulation time. Recall that R is to be distinguished from Rg(t), Rg, RGg , and REg defined earlier. The running time average of R, Rg(t), is reported in panel b. Its long time value Rg can be read from the long time terminus of the curve in the middle panel. A Kratky plot of the scattered intensity is shown in panel c. The points are the experimental data,6 and the dashed curve was computed from the long time running average of P(q) as calculated using I(q) in eq 1. The numerical values of REg , Rg, and RGg are given in Table 1. To permit easy comparison with the earlier RIS results,6 RRIS g is reported in column 5 of Table 1. The latter quantity was computed from the RIS calculation using the Guinier

approximation as described above for computing RGg . Values of REg reported in Table 1 differ slightly from those reported earlier6 because we chose here to use a consistent procedure for extracting the radius of gyration from the experimental and theoretical results using the Guinier approximation. Experimental uncertainties reported for REg are all about (2.5%. To convey the range of radii of gyration characteristic of an ensemble of G3 molecules, column 6 of Table 1 reports the standard deviation, σR, of the distribution of values of R shown in panel a. For G3, Rg and RGg are not significantly different, so either is appropriate for comparison with REg . For the larger oligomers Rg and RGg become increasingly different (Table 1), suggesting that the experimental range of q did not descend adequately into the limiting range of q for which the Guinier approximation is valid. Hence, we will focus primarily on making comparisons among REg , RGg , and RRIS g , all of which are based on using the Guinier approximation. The trisaccharide G3 contains two (1f4)-linkages and no (1f6)-linkages. Thus, molecular flexibility is governed, ignoring any possible interactions of greater range than nearest neighbor residues, by the maltose conformational energy surface in Figure 4. The range of fluctuations in R is about (5% and the standard deviation σR (Table 1) is about (2%. Thus, G3 is a relatively rigid molecule, the (1f4)linkages permit little latitude for conformational fluctuations and the six-membered sugar rings themselves are quite rigid, so the earlier RIS approximation of a single rotational isomeric state for the (1f4)-linkage seems again quite well justified. The MD simulation fits REg well within experi-

1244

Biomacromolecules, Vol. 6, No. 3, 2005

Jaud et al.

Figure 7. Results of the AMBER*+GB/SA simulation for G3. (a) Instantaneous radius of gyration R, (b) running time average rootmean-square radius of gyration Rg(t), (c) Kratky plots: the dots represent experimental results, whereas the plain line is the result of the simulations. Table 1. Rg Results for AMBER* Simulations of Pullulana

Rg (Å)

SAXS sample

REg

Rg

RG g

RRIS g

〈R〉 ( σR

G3 (G3)2 (G3)3 (G3)4

4.56 ( 0.11 6.56 ( 0.16 8.03 ( 0.20 9.88 ( 0.25

4.62 ( 0.01 6.63 ( 0.08 8.74 ( 0.08 11.38 ( 0.10

4.60 6.50 8.16 10.78

4.56 6.61 8.17 9.87

4.62 ( 0.09 6.60 ( 0.64 8.70 ( 0.77 11.36 ( 1.04

a RE is the experimental R and its uncertainty; R refers to R (∞) g g g g calculated from the atomic coordinates (eq 3) and its uncertainty is the G RIS standard deviation of Rg(t); Rg and Rg refer to Rg calculated using the MD and RIS models (and the Guinier approximation). The last column shows the average value of R, the instantaneous radius of gyration, and its standard deviation σR.

mental error, so the model accurately describes the global dimensions of the pullulan repeating unit corresponding to the region of small q. This is, however, not a very demanding test of the model, and we have easily fit REg equally well with other approaches.6,10 The fit to the scattering function at higher q is also quite adequate, although perhaps not quite so good as achieved with the RIS approach.6,10 Note that the scatter of the Kratky data at the upper end of the experimental q range arises from the smaller signal-to-noise ratio observed at large angle. There is also a greater likelihood of systematic errors in the high angle data due to corrections and calibrations that must be applied.6

Figure 8. Results of the AMBER* + GB/SA simulation for (G3)2. (a) Instantaneous radius of gyration R, (b) running time average rootmean-square radius of gyration Rg(t), (c) Kratky plots: the dots represent experimental results, whereas the plain line is the result of the simulations.

Pullulan Oligomer (G3)2. Figure 8 shows a 100 ns simulation of the hexameric pullulan oligomer (G3)2. The results are similar to those reported for G3 in Figure 7, but the amplitude of fluctuations in R is much greater for (G3)2 (panel a and Table 1) and R does not fluctuate around a single value as it does for G3. The latter feature of the R curve gives rise to two segments on the curve for Rg(t) (panel b) with positive slope. Between 41 and 50 ns and between 98 and 99 ns the hexamer adopts a conformation distinctly more extended than the one that dominates during most of the simulation. The first departure from the more compact conformation, because it comes relatively early in the simulation and has a significant duration, has a larger impact on the running average Rg(t) than does the later departure. The relatively infrequent occurrence of the more extended conformer during the simulation naturally raises the issue of whether the long time running average Rg reflects ergodic sampling. This can only be answered definitively by extending the simulation significantly. Nevertheless, the simulation results for (G3)2 fit consistently with those for the other three oligomers, and we proceed under the assumption that Rg for (G3)2 is not seriously compromised by poor sampling statistics. The greater amplitude of fluctuations in R for (G3)2 compared to G3 is expected simply on the basis of the larger

Biomacromolecules, Vol. 6, No. 3, 2005 1245

Simulation of Pullulan Oligomers

Table 3. List of the R Values Adopted by the Pullulan Hexamer (G3)2 while in the Minimum of Each of the Wells on the Potential Surface of Isomaltose (Figure 5)a (G3)6 conformation

R (Å)

A4-A6 A4-B6 A4-C6 A4-D6 A4-E6 A4-F6

7.08 6.78 4.67 7.68 7.46 6.77

a The (φ ,ψ ) angles were set to correspond to minimum A on the 4 4 4 potential surface of maltose (Figure 4), and φ6 was set to 60° to meet the requirements of the exo-anomeric effect.

Figure 9. Conformations that the pullulan hexamer adopts at certain values of its radius of gyration. Only carbon and oxygen atoms are shown. Table 2. How the Glycosidic Linkage Dihedrals of the Pullulan Hexamer Change with the Radius of Gyrationa linkage

dihedral

R ∼ 5.8

R ∼ 6.5

R ∼ 7.4

R ∼ 8.0

1 (1f4) 2 (1f4) 3 (1f6)

φ41 ψ41 φ42 ψ42 φ63 ψ63 ω63 φ44 ψ44 φ45 ψ45

84 -143 95 -153 53 -133 76 86 -130 95 -138

77 -154 61 -157 58 -158 59 74 -160 87 -142

70 173 98 -147 58 175 30 89 -137 71 -171

97 -147 81 -145 63 176 -62 60 -145 93 -156

4 (1f4) 5 (1f4)

a The first linkage is the (1f4)-linkage uniting the terminal sugar at the nonreducing end to the adjacent sugar ring; the third linkage in (G3)2 is the sole (1f6)-linkage in the oligomer; the fifth linkage is the (1f4)linkage connecting the terminal (reducing) ring to the preceding ring. The terminology used is better explained with an example: φ41 designates the φ4 dihedral located in linkage 1.

overall dimensions of the higher oligomer. This effect is accentuated by the occurrence of the (1f6)-linkage in the center of the hexamer (Figure 1). Figure 9 shows four conformations of (G3)2 extracted from the simulation with their respective values of the instantaneous radius of gyration R. The most extended and most compact of these have values of R close to the extrema displayed in Figure 8a. The conformer with R ) 6.5 Å has dimensions very close to the long time average Rg ) 6.63 Å (Table 1). Table 2 reports the values of the glycosidic dihedral angles in the four conformers of (G3)2 shown in Figure 9. These angles are designated, for example, by ψ42, where the first subscript designates the type of linkage and the second identifies the position of that linkage in the sequence of five linkages in

the hexamer, counting from the nonreducing end. All five of the angles φij are in the narrow range 50-100° dictated by the exo-anomeric effect. The angles ψ4j are likewise confined to the 55° range from -130 to +175° as expected from the conformational energy map in Figure 4. For none of these angles (φij, ψ4j) is any systematic correlation with the value of R evident in Table 2. There is, however, a strong correlation of R with the value of ψ63, which adopts a conformation that is increasingly close to the 180° trans state as R increases. At the same time, ω63 changes monotonically from +76° to -62° as R becomes larger. It is clear that the most extended conformations of the hexamer is approximately (ψ63, ω63) ) (180°, -60°), which corresponds to the highly populated minimum F6 in Figure 5 for isomaltose. Table 3, however, reveals that for characteristic values of (φ4,ψ4, φ6), the conformations included in this minimum can have length comparable to the ones in B6. In Figure 9, the more compact forms of the hexamer are characterized by (ψ63, ω63) conformations corresponding to well-populated minimum B6 in Figure 5 and, especially, the domain between minima B6 and C6. It is noteworthy that Liu et al. concluded that substantial occupation of minimum C6, which produces very compact conformations of (G3)2 (R ) 4.67 Å), was essential for consistency of their RIS model for pullulan with experimental SAXS data.6 Given the relative inflexibility of the (1f4)-linkages, it is no surprise that the range of values of R reported in Figures 8 and 9 is due almost entirely to differences in the torsion angles associated with the central (1f6)-linkage, in particular differences in ω6 and ψ6, since φ6 is rather tightly constrained by the exo-anomeric effect. Values of Rg and RGg from the (G3)2 simulation are again in excellent agreement with the experimental result (Table 1); that is, the fit to the Kratky plot at low q is essentially perfect. The fit at high q falls below the experimental points, as was the case for the G3 simulation (Figure 7c). An additional factor enters the effort to match the experimental results for (G3)2. Because the experimental sample of (G3)2 was contaminated with a small amount of G3, it is necessary to make an appropriately weighted average of the results from the simulations for (G3)2 and G3 in order to make the comparison with the experimental results. All the results from the simulation displayed in Figure 8 include a correction for the contribution from the G3 contaminant. Use of the z-distribution to correct the theoretical values of P(q) and the associated radii of gyration is straightforward and is explained elsewhere.6 A similar procedure is required in order to compare the results of the simulations of (G3)3 and (G3)4

1246

Biomacromolecules, Vol. 6, No. 3, 2005

Jaud et al.

Figure 10. Results of the AMBER* + GB/SA simulation for (G3)3. (a) Radius of gyration, (b) root-mean-square radius of gyration, (c) Kratky plots: the dots represent experimental results, whereas the plain line is the result of the simulations.

Figure 11. Results of the AMBER* + GB/SA simulation for (G3)4. (a) Radius of gyration, (b) root-mean-square radius of gyration, (c) Kratky plots: the dots represent experimental results, whereas the plain line is the result of the simulations.

with experiment, because both of these oligomers were contaminated with each of the lower oligomers. Pullulan Oligomer (G3)3. Results of the simulation of (G3)3 in Figure 10 show that R departs several times during the 100 ns simulation from its dominant value. The maximum excursion of the fluctuations in R is greater than for (G3)2, which is consistent with the greater length of the oligomer backbone. The running average R(t) rises in several steps as the simulation progresses. It is not obvious that Rg has truly achieved constancy even after 100 ns, but the rate of change with time is clearly approaching zero. The fit to the Kratky plot of the experimental data appears to be better than for either G3 or (G3)2, This is certainly true for the regime of large q, but examination of Table 1 shows that Rg exceeds REg substantially. On the other hand, in this case Rg > RGg and RGg matches REg , also extracted using the Guinier approximation, within experimental error. Likewise, RRIS g also matches REg just as well. Indeed, the MD simulations and the RIS model, both based on the same AMBER* + GB/SA force field, are equally competent at matching RGg with REg for the oligomers G3, (G3)2, and (G3)3, while the RIS approach is perhaps marginally better at fitting the higher angle scattering data. Pullulan Oligomer (G3)4. Plots of simulation results for R and Rg(t) in panels a and b of Figure 11 share the characteristics of their counterparts for the smaller oligomers.

It again appears that Rg corresponding to 100 ns has achieved a value close to the true infinite time result. However, as the chain length of the oligomer increases, with a concomitant increase in the total number of accessible conformations, it is less easy to be persuaded that 100 ns is long enough to sample all of the configuration space ergodically, so the simulation was extended to 200 ns. Panel b, however, indicates that Rg does not vary significantly after 100 ns, suggesting that extending the simulation was not absolutely required. The original G3, (G3)2, and (G3)3 simulations, which we need to incorporate to take the shorter oligomers contamination into account, were 100 ns too short to be used with the 200 ns (G3)4 run. Rather than extending them, we made the reasonable assumption that they reached ergodic sampling and duplicated them once so their total length was 200 ns. These extended simulations were then used to generate the second half (t > 100 ns) of Figure 11. The fit to the experimental Kratky plot in panel c of Figure 11 is perceptibly worse than those achieved with the other three lower oligomers. Negative departure of the simulated result from the experimental data at relatively small q signifies that the simulation will substantially overestimate REg (Table 1), although the simulation does catch up with the experimental data as q approaches 0.5. Here Rg exceeds REg by 15% and even RGg is too large by 9%. This is to be E contrasted with the RIS result RRIS g , which matches Rg

Simulation of Pullulan Oligomers

Figure 12. Trajectories for the aqueous pullulan hexamer (G3)2 (AMBER*). Contour lines represent increments of 1 kcal/mol.

exactly. Thus, two important observations emerge. First, RGg departs increasingly from Rg as the chain length of the oligomer increases, presumably because the Guinier approximation is increasingly inaccurate with increasing oligomer size. Second, the MD simulation for (G3)4 predicts a substantially larger value of the radius of gyration than is observed and then is predicted by the RIS treatment. This is true even when the comparison is made between RGg and REg , which are equally dependent on the Guinier approximation. Evidence of Excluded Volume Effects in Higher Pullulan Oligomers. The observation that RGg > RRIS g is almost certainly a reflection of the influence of the excluded volume effect, which must come increasingly into play in the MD simulation as the oligomer chain length grows. One can get some insight into this issue by examining the distribution among the states of the linkage torsion angles for oligomers of different sizes. Figure 12 shows the distribution among the torsional states of the single (1f6)-linkage in the pullulan hexamer (G3)2 during the 100 ns simulation reported in

Biomacromolecules, Vol. 6, No. 3, 2005 1247

Figure 8. It is most instructive to compare Figure 12 with the distribution among the torsional states of the (1f6)linkage in isomaltose reported in Figure 5. To facilitate this comparison with the isomaltose simulation, which lasted only 48 ns, we have sampled the trajectories at the same number of points (9600) in Figures 5 and 12. It is evident that minimum C6, which leads to very compact conformations of the hexamer, is never visited in the hexamer simulation, whereas it shows reasonable population in isomaltose. At the same time, that portion of minimum A6 in the range 30 e ψ63 e 60° is less heavily populated in the hexamer than in isomaltose. There is also transfer of density out of minima E6 and F6 for the hexamer in comparison to isomaltose. The density missing in Figure 12 from minima C6, A6, E6, and F6 appears instead in the very extended minimum D6, which is much more heavily populated in the hexamer than in isomaltose, and, especially, in minimum B6. Plots in panels c and d of Figure 13, made respectively for the sixth and ninth (1f6)-linkages in (G3)4, are essentially indistinguishable from panel a of Figure 12. The same plot in panel b of Figure 13 for the third (1f6)-linkage of (G3)4, i.e., that closest to the nonreducing end of the chain, shows somewhat less transfer of density to D6 than in panel a of Figure 12, but the differences are slight. The data in Figures 12 and 13 are presented numerically in Table 4, where it is seen that in absolute terms (1f6)-linkage population density in the oligomers is transferred most heavily to minimum B6 in comparison to isomaltose, while in relative terms minimum D6 receives the largest population increment. Figure 14 shows the torsional angle distribution for (1f4)linkages 2, 5, 8, and 11 in (G3)4. Comparison with Figure 4 for the (1f4)-linkage of maltose shows that incorporation of this linkage into pullulan has very little effect on its torsional state distribution. Minimum B4 was visited occasionally by (1f4)-linkage 2, as was the case for maltose, but this subsidiary minimum was avoided altogether by (1f4)-linkages 5, 8, and 11. Linkage 11 visits the region ψ411 > 270° more often than do the other (1f4)-linkages. Linkages 2 and 11 are near the chain termini, and somewhat greater conformational latitude is to be expected. Figure 15 shows the pullulan nonamer, (G3)3, with successive (1f6)-linkages in conformation C6. The ball-andstick model (Figure 15a) is projected to minimize the impression of steric encounters between components of the chain. When realistic van der Waals radii are added to the atoms (Figure 15b), the highly congested nature of this conformation is evident. Minimum C6 is avoided even in the hexamer, so this allowed conformation of isomaltose is already problematic in oligomers smaller than (G3)3. Potential difficulties with C6 were noted in the RIS treatment of pullulan, where it was observed that adoption of the C6 conformation for the (1f6)-linkage of the pullulan fragment panose, O-R-D-glucopyranosyl-(1f6)-O-R-D-glucopyranosyl-(1f4)-D-glucose, provokes some close contacts between second neighboring glucose residues in the panose fragment.6 It was nevertheless necessary to include minimum C6 with a probability of 0.158 in order to fit the SAXS data on pullulan oligomers and the light scattering data on the high molecular weight polymer.

1248

Biomacromolecules, Vol. 6, No. 3, 2005

Jaud et al.

Figure 13. Trajectories for (a) isomaltose and for the (b) third, (c) sixth, and (d) ninth (1f6)-linkages (starting from the nonreducing ring extremity) of the aqueous pullulan dodecamer (G3)4 (AMBER*). Contour lines represent increments of 1 kcal/mol. Table 4. Percentage of Occupation of the Different Minima of the (1f6)-Linkages of Isomaltose and Pullulan Oligomers location

percentage of occupation

minimum

[ψ6]/[ω6]

isomaltose

(G3)2

(G3)4 3

(G3)4 6

(G3)4 9

A6 B6 C6 D6 E6 F6

[0°-120°]/[0°-120°] [120°-260°]/[0°-120°] [260°-340°]/[40°-120°] [90°-230°]/[120°-220°] [40°-130°]/[-110°-0°] [130°-250°]/[-110°-0°]

31.1 39.4 4.0 1.9 6.9 16.4

4.5 70.4 0.0 12.8 3.1 8.8

9.1 62.6 0.0 6.4 5.7 16

7.7 54.3 0.0 18.5 4.9 14.4

7.0 58.3 0.0 15.9 4.9 13.6

a Location describes the (ψ ,ω ) domain associated with each minimum. (G ) corresponds to the pullulan hexamer, whereas (G ) n refers to the 6 6 3 2 3 4 (1f6)-linkage that is the nth linkage of pullulan dodecamer counting from the nonreducing end.

Even in isomaltose, however, minimum C6 makes a relatively small contribution to the mean spatial properties of the molecule. Understanding the progression of (1f6)linkage population distributions with increasing oligomer length requires, more importantly, understanding the depletion of minimum A6 from 31.1% occupancy in isomaltose to something between 4.5 and 9.1% occupancy in the pullulan oligomers, and the less dramatic depletion of density in states E6 and F6 (Table 4). Figure 16 shows that occupancy of minimum A6 in the hexamer (G3)2 provokes a short contact in the panose fragment between the pendant hydroxymethylene group of the reducing residue of that fragment and the nonreducing sugar. This short contact between sugars 3 and 5 of the hexamer does not arise in the disaccharide isomaltose nor in the trimer G3, which contains no (1f6)linkages. Consequently, the population of A6 is severely

depleted in the pullulan oligomers beyond G3. It can be noted, however, from Figures 5 and 12 that only the portion of minimum A6 at the lower values of ψ6i is excluded by the close contact in question. The same behavior is observed for (G3)4 in Figure 13. Much of the density that shifts out of A6 ends up in B6, which, as long as the region where it merges with C6 is avoided, corresponds to moderate extension of the (G3)2 hexamer with radii of gyration close to the ensemble mean value Rg (Tables 1 and 3) and freedom from any unacceptably close contacts. Population density transferred out of C6 and A6 in isomaltose goes secondarily to minimum D6. Examination of Figures 2 and 5 shows that minima B6 and D6 differ by rotation of ω6 through approximately 120° about the C5-C6 bond in the pullulan backbone. The chain conformations of the hexamer (G3)2 in both of these minima are similar, are

Simulation of Pullulan Oligomers

Biomacromolecules, Vol. 6, No. 3, 2005 1249

Figure 15. Pullulan nonamer in the C6 conformation. Figure 14. Trajectories for the (a) second, (b) fifth, (c) eighth, and (d) eleventh (1f4)-linkages of aqueous pullulan dodecamer (G3)4 (AMBER*). Contour lines represent increments of 1 kcal/mol.

free of high energy self-intersections of the chain, and have radii of gyration near Rg or slightly larger (Tables 1 and 4). Hence, it is clear why minimum D6 serves as the secondary repository of population density transferred out of C6 and A6 as one moves from isomaltose to pullulan oligomers containing (1f6)-linkages. It is not so clear why there is some shift of density out of minima E6 and F6 on moving from isomaltose to pullulan. Figures 12 and 13 clearly show that the contour line encompassing most of the trajectories of these minima is one level higher for isomaltose than for the hexamer (G3)2, but (G3)2 in these minima still has an extension comparable in length to B6 and D6, and there are no obvious close contacts. because the RIS treatment So it appears that RGg > RRIS g ignores interactions beyond first neighbor sugar residues and therefore includes a finite likelihood of conformations that will be self-intersecting at the range of a few sugars. If the criterion for success of the model is simply that it fits the SAXS data, the nearest neighbor RIS treatment6 works better than the presumably more realistic MD approach. This conclusion was also reached using an extended RIS approach that accounts for longer range interactions in the pullulan oligomers.10 The present simulations show, however, that the success of the RIS approach must be based in part on compensating errors in the treatment.

Figure 16. Pullulan hexamer with central (1f6)-linkage in the A6 conformation. All (1f4)-linkages are in conformation A4. The panose fragment involves residues 3, 4, and 5 counting from the left (nonreducing) end.

Conclusion Because the RIS6 and MD treatments of the pullulan oligomers were both carried out with AMBER* and GB/SA solvation, we have here an opportunity to compare these two methods directly under the influence of a common force field. The RIS approach6 is computationally efficient and lends itself to ready physical interpretation. Moreover, it appears to provide a good fit to experimental scattering data for all four pullulan oligomers and for the pullulan high polymer. Although the MD method yields agreement between Rg and REg within experimental uncertainty for G3, (G3)2, and (G3)3,

1250

Biomacromolecules, Vol. 6, No. 3, 2005

this is not the case for (G3)4, which is found in the simulation to be too extended relative to the experimental result. This discrepancy is absent in the RIS results because that treatment allowed (indeed, required) a nonnegligible contribution from compact (1f6)-linkage conformations (corresponding to to fit REg . We have minimum C6) in order to force RRIS g shown here that, in the course of adjusting RIS model parameters to obtain the best fit to the experimental SAXS observations on the pullulan oligomers,6 features were introduced into the model that allowed for local conformations that are physically unrealistic due to excluded volume effects at a range in the chain sequence greater than that encompassed by the RIS model, which is built on dimeric segments. These excluded volume effects, which come into play initially at the level of the hexamer or, indeed, at the level of the trimeric segment panose,6 are clearly disclosed here using the more computationally intensive MD approach. Although the more realistic MD approach is less computationally efficient than the RIS approach, we have demonstrated here that useful results can be obtained without inclusion of an explicit solvation model. The 100 ns simulation for the dodecamer (G3)4 required approximately 1 month with the relatively modest hardware configuration described above. On the other hand, we estimate from results not reported here in detail that a 100 ns simulation for the hexamer (G3)2 using an all-atom force field and the same hardware would require 1 order of magnitude more computing time. Obviously, the high speed provided by using the GB/SA or other continuum solvation model can provide a considerable advantage in achieving ergodic sampling of the system. The possible sacrifice in accuracy associated with the use of a continuum solvation model might be more than compensated by the shorter times needed to achieve ergodic sampling in simulations designed, for example, to test the ability of the polymer chain model to reproduce nano- or picosecond molecular dynamics data (e.g., from NMR relaxation or neutron scattering experiments).39-45 We note finally that data on frequency of rotational isomeric state occupation of the sort presented in Figures 13 and 14 and Table 4 can provide the basis for constructing mean-field potentials for use in the spirit of the “nearly separable configuration energy RIS approach” of Stokke et al.46 Although we have disclosed here the flaws in a computationally efficient RIS approach to the modeling of polymer conformational statistics based on relatively short polymer chain segments, we have not explained to this point why the current MD simulation has failed to reproduce the experimental REg values for the longest of the oligomers studied, (G3)4. This is presumably related to inaccuracies in the force field employed. Acknowledgment. This work was funded by NIH grant GM33062 to D.A.B. and by NSF grant MCB-0078278 to D.J.T. The authors thank the Quebec “Fonds pour les Chercheurs et l′Aide a` la Recherche” for a graduate fellowship for S.J.

Jaud et al.

References and Notes (1) Yalpani, M., Polysaccharides: syntheses, modifications, and structure/ property relations, 1st ed.; Elsevier: Amsterdam, 1988. (2) Tsujisaka, Y.; Mitsuhashi, M. Pullulan. In Industrial Gums: Polysaccharides and Their DeriVatiVes, 3rd ed.; Whistler, R. L., BeMiller, J. N., Eds.; Academic Press: San Diego, CA, 1993; pp 447460. (3) Morris, V. J. In Food Polysaccharides and Their Applications, 1st ed.; Stephen, A. M., Ed.; Marcel Dekker: New York, 1995; pp 341375. (4) Colson, P.; Jennings, H. J.; Smith, I. C. P. J. Am. Chem. Soc. 1974, 96, 8081-8087. (5) Buliga, G. S.; Brant, D. A. Int. J. Biol. Macromol. 1987, 9, 7176. (6) Liu, J. H.-Y.; Brant, D. A.; Kitamura, S.; Kajiwara, K.; Mimura, M. Macromolecules 1999, 32, 8611-8620. (7) Brant, D. A.; Burton, B. A. In Solution Properties of Polysaccharides, 1st ed.; Brant, D. A., Eds.; American Chemical Society: Washington, DC, 1981; Vol. 150, pp 81-99. (8) Burton, B. A.; Brant, D. A. Biopolymers 1983, 22, 1769-1792. (9) Buliga, G. S.; Brant, D. A. Int. J. Biol. Macromol. 1987, 9, 7786. (10) Liu, J. H.-Y.; Brameld, K. A.; Brant, D. A.; Goddard, W. A., III. Polymer 2000, 43, 509-516. (11) Qiu, D.; Shenkin, P. S.; Hollonger, F. P.; Still, W. C. J. Phys. Chem. A 1997, 101, 3005-3014. (12) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. J. Am. Chem. Soc. 1990, 112, 6127-6129. (13) Mohamadi, F.; Richards, N. G. J.; Guida, W. C.; Liskamp, R.; Lipton, M.; Caulfield, C.; Chang, G.; Hendrickson, T.; Still, W. C. MacroModel-An Integrated Software System for Modeling Organic and Bioorganic Molecules Using Molecular Mechanics. J. Comput. Chem. 1990, 11, 440-467. (14) Honig, B.; Nicholls, A. Science 1995, 268, 1144-1149. (15) Tannor, D. J.; Marten, B.; Murphy, R.; Friesner, R. A.; Sitkoff, D.; Nicholls, A.; Ringnalda, M.; Goddard, W. A., III; Honig, B. J. Am. Chem. Soc. 1994, 116, 11875-11882. (16) Ringnalda, M. N.; Langlois, J.-M.; Greeley, B. H.; Murphy, R. B.; Russo, T. V.; Cortis, C.; Muller, R. P.; Marten, B.; Donnelly, R. E.; Mainz, D. T.; Wright, J. R.; Pollard, W. T.; Cao, Y.; Won, Y.; Miller, G. H.; Goddard, W. A., III; Friesner, R. A. Jaguar 3.0; Schro¨dinger, Inc.: Portland, OR, 1997. (17) Tobias, D. J. Curr. Opin. Struct. Biol. 2001, 11, 253-261. (18) Ha, S. N.; Giammona, A.; Field, M.; Brady, J. W. Carbohydr. Res. 1988, 180, 207-221. (19) Brady, J. W.; Schmidt, R. K. J. Phys. Chem. 1993, 97, 958966. (20) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey;, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (21) Stoddart, J. F., Stereochemistry of Carbohydrates, 1st ed.; WileyInterscience: New York, 1971. (22) Recommendations (1), Conformational nomenclature for five- and six-membered ring forms of monosaccharides and their derivatives. Eur. J. Biochem. 1980, 111, 295-298. (23) Recommendations (4), Abbreviated terminology of oligosaccharide chains. Eur. J. Biochem. 1982, 126, 433-437. (24) Recommendations (5), Polysaccharide nomenclature. Eur. J. Biochem. 1982, 126, 439-441. (25) Recommendations (6), Symbols for specifying the conformation of polysaccharides. Eur. J. Biochem. 1983, 131, 5-7. (26) Pe´rez, S.; Marchessault, R. H. Biopolymers 1979, 18, 23692374. (27) Nishida, Y.; Hori, H.; Ohrui, H.; Meguro, H. Carbohydr. Chem. 1988, 7, 239-250. (28) Bock, K.; Duus, J. O. J. Carbohydr. Chem. 1994, 13, 513-543. (29) Ohrui, H.; Nishida, Y.; Watanabe, M.; Hori, H.; Meguro, H. Tetrahedron Lett. 1985, 26, 3251-3254. (30) Wolfe, S. Acc. Chem. Res. 1972, 5, 102-111. (31) Cramer, C. J.; Truhlar, D. G.; French, A. D. Carbohydr. Res. 1997, 298, 1-14. (32) Ha, S. N.; Madsen, L. J.; Brady, J. W. Biopolymers 1988, 27, 19271952. (33) Stortz, C. A. Carbohydr. Res. 1999, 322, 77-86. (34) Tobias, D. J., In Hydration Processes in Biology: Theoretical and Experimental Approaches, 1st ed.; Bellissent-Funel, M.-C., Eds.; IOS Press: Amsterdam, 1999; pp 379-409. (35) Andersen, H. C. J. Comput. Phys. 1983, 52, 24-34.

Simulation of Pullulan Oligomers (36) Glatter, O.; Kratky, O. Small Angle X-ray Scattering, 1st ed.; Academic Press: London, 1982. (37) Wilson, A. J. C., International Tables for Crystallography; The International Union of Crystallography; Kluwer Academic Publishers: Boston, 1992; C: Mathematical, Physical and Chemical Tables. (38) Guinier, A. X-ray Diffraction in Crystals, Imperfect Crystals, and Amorphous Bodies, 1st ed.; W. H. Freeman and Company: San Francisco, CA, 1963. (39) Brant, D. A.; Liu, H.-S.; Zhu, Z. S. Carbohydr. Res. 1995, 278, 1126. (40) Brant, D. A. Pure Appl. Chem. 1997, 69, 1885-1892. (41) Perico, A.; Mormino, M.; Urbani, R.; Cesa`ro, A.; Tylianakis, E.; Dais, P.; Brant, D. A. J. Phys. Chem. B 1999, 103, 8162-8171.

Biomacromolecules, Vol. 6, No. 3, 2005 1251 (42) Cavalieri, F.; Chiessi, E.; Paci, M.; Paradossi, G.; Flaibani, A.; Cesa`ro, A. Macromolecules 2001, 34, 99-109. (43) Letardi, S.; La Penna, G.; Chiessi, E.; Perico, A.; Cesa`ro, A. Macromolecules 2002, 35, 286-300. (44) Corzana, F.; Motawia, M. S.; Herve´ du Penhoat, C.; Perez, S.; Tschampel, S. M.; Woods, R. J.; Engelsen, S. B. J. Comput. Chem. 2004, 25, 573-586. (45) Smith, L. J.; Price, D. L.; Chowdhuri, Z.; Brady, J. W.; Saboungi, M.-L. J. Chem. Phys. 2004, 120, 3527-3530. (46) Stokke, B. T.; Talashek, T. A.; Brant, D. A. Macromolecules 1994, 27, 1124-1135.

BM049463D