Analysis of Solvation and Gelation Behavior of Methylcellulose Using

Article pubs.acs.org/JPCB

Analysis of Solvation and Gelation Behavior of Methylcellulose Using Atomistic Molecular Dynamics Simulations Wenjun Huang, Indranil S. Dalal,§ and Ronald G. Larson* Department of Chemical Engineering, University of Michigan, Ann Arbor, Michigan 48109-2136, United States S Supporting Information *

ABSTRACT: We adopt atomistic molecular dynamics (MD) simulations to study the solvation and gelation behavior of homogeneous methylcellulose (MC) and model random oligomers that represent the commercial cellulosic polymer product METHOCEL A in water and acetone solvents. We demonstrate that the two carbohydrate-specific GROMOS force fields, GROMOS 45A4 and 56Acarbo, are capable of reproducing characteristic angle distributions and the persistence length of MC chains reported in the literature. We then use the GROMOS 56Acarbo force field in both single-chain and multiple-chain simulations to characterize their solvation behavior through radial distribution functions, hydrogenbond counts, and contact map analyses. We find that the un-methylated O6 position on the cellulose ring forms the most hydrogen bonds, followed by O2 and O3, implying that methylation at the 6 position reduces hydrogen bonding more than does methylation at other positions. O6−O6 is the most probable intermolecular hydrogen bond between different MC molecules. Dimethylated and trimethylated MCs form aggregated structures at low temperatures and precipitate-like structures at high temperatures in water but disperse randomly in acetone. This is consistent with experimental observations of gelation at elevated temperatures in water. The heterogeneous METHOCEL A model shows increased aggregation of trimethylated monomer units at elevated temperatures, suggesting that hydrophobic interaction is the main factor that induces the gelation, rather than hydrogen bonding.

■

INTRODUCTION Methylcellulose (MC) and its derivatives, including hydroxypropyl methylcellulose, are categorized as safe to use as food additives by the U.S. Food and Drug Administration.1 The METHOCEL products from the Dow Chemical Company, which are based on methylcellulose, are also widely used as adhesives, agricultural chemicals, ceramics processing aids, construction products, and pharmaceutical applications.2 Besides its great commercial value, methylcellulose is of considerable research interest both experimentally and theoretically.3,4 The monomer unit of cellulose is the anhydroglucose unit (AGU) (see Figure 1 for details). When one or more hydroxyl groups (−OH) on AGUs are substituted with methyl groups (−CH3), the product is methylcellulose, and AGUs become methylcellulose monomers. Because each methylcellulose monomer has three functional groups available for substitution, methylcellulose can be made with a range of hydrophobicities by varying the degree of substitution (DS), defined as moles of substituents (i.e., hydrophobic methyl groups) per mole of methylcellulose monomer.5 Although each AGU has three hydroxyl groups, natural cellulose is insoluble in water because the hydroxyl groups form an extensive intermolecular hydrogen-bonding network.6 Methylcellulose with fully substituted hydrophobic functional groups (i.e., methylcellulose with DS = 3), on the other hand, is also insoluble in water because of its strong hydrophobic interactions. Interestingly, partially substituted hydrophobically modified methylcellulose (Figure 1) with DS ≈ 2 is water© XXXX American Chemical Society

soluble and is therefore of great commercial and scientific value.5 Although the structural properties of these methylcelluloses in water and other solvents have been extensively characterized over the past few decades,5,7−10 a complete picture of their solvation behavior has not been revealed. One of the important properties of methylcellulose is that it forms a thermoreversible gel upon heating. Although this transition has been well-studied, the actual gel-forming mechanism is still not clearly understood. A key question that any mechanism must address is how intermolecular hydrophobic interactions form a network and avoid simple phase separation. Although several competing mechanisms have been proposed, all involve nonuniformity of hydrophobic substitution in some way. Kato et al.11 proposed that trisubstituted methylcellulose units, which occur as a result of the randomness of substitution, act as cross-link junctions. Li et al.12 suggested that the more hydrophobic regions, again due to heterogeneities in substitution, associate at elevated temperature to form cross-linked structures. Haque and Morris,13 on the other hand, proposed that methylcellulose chains in solution are bundled with parallel chain configurations at low temperature. Upon heating, the unsubstituted regions remain intact because of the strong hydrogen bonding present in natural cellulose while the hydrophobically modified regions Received: September 26, 2014 Revised: November 6, 2014

A

dx.doi.org/10.1021/jp509760x | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

the formation of a thermoreversible gel network using atomistic molecular dynamics (MD) simulations. To our knowledge, our study is the first to examine intra- and intermolecular structure formation, hydrogen bonding, and clustering in solutions of both homogeneously and heterogeneously substituted methylcellulose oligomers using atomistic MD simulations. The model heterogeneously substituted oligomers we develop are based on the substitution statistics of METHOCEL A, a commercial methylcellulose product from the Dow Chemical Company, so our results are of practical interest. The paper is organized as follows. We first discuss and validate our choice of force fields. We then discuss chain flexibility and hydrogen-bond statistics for a single cellulosic chain in water. Next we discuss hydrogen bonding, molecular radial distribution functions, and oligomer clustering of multiple homo-oligomers and model heterooligomers based on the METHOCEL A chemistry. Throughout we emphasize the effect of temperature, since it plays a crucial role in the gelation of methylcellulose.

Figure 1. Schematic of cellulose and methylcellulose. The anhydroglucose unit (AGU) is the monomer unit of cellulose, where each R group is a hydrogen. In monomers of methylcellulose, one or more R groups are methyl groups (−CH3). Both AGUs and methylcellulose monomers are connected via β(1→4) linkages. The primes indicate atoms in the second (right) monomer in any two consecutive monomers in the cellulosic chain, although the two monomers are bonded identically if they are not terminal monomers (see the text for details). O2, O3, and O6 are the three reactive hydroxyl groups (−OH) on each AGU that can be substituted with different functional groups (i.e., methyl groups) so that the AGU becomes a methylcellulose monomer. The blue dashed lines indicate the two most predominant intramolecular hydrogen bonds that are present in cellulose.

■

METHODS Methylcellulose monomers with all possible combinations of methyl substituents were built using Materials Studio, version 6.1 (Accelrys Inc.), which is a commercial licensed software package. These monomers were then used to build homooligomers and random oligomers with user-defined probabilities for incorporation of each monomer. Two different force fields were evaluated, namely, GROMOS 56Acarbo,26 which is a more recent force field, and the older GROMOS 45A4,27 from which the GROMOS 56Acarbo is derived. Both force fields are carbohydrate-specific, with the newer one featuring an improved treatment of the ring atoms in the cellulosic repeat unit. The simulations were carried out in GROMACS, version 4.5.5.28−30 For the single-chain simulations, cubic boxes with side lengths comparable to the end-to-end distance of the oligomer chains were constructed. The box sizes for the singlechain simulations ranged from 6 to 22 nm. The simulations of multiple homo-oligomers were carried out in a 12 nm cubic periodic box. The simulations based on the METHOCEL A chemistry were carried out in cubic boxes with sizes between 20 and 22 nm. The concentrations of the oligomers ranged from 2 to 6 wt %. Unless otherwise stated, the oligomers were solvated with simple point charge (SPC) water.31 The density of the initial simulation box was approximately equal to the bulk water density (1 g/cm3). All of the systems were subjected to 20 000 steps of energy minimization using a steepest-descent method. A 0.5 ns NVT equilibration followed by a 10 ns NPT equilibration was conducted for each system using a time step of 1 fs. A weak temperature coupling using a velocity-rescaling thermostat32 with a relaxation time of 0.2 ps and a weak pressure coupling using a Berendsen barostat33 with a relaxation time of 0.5 ps were used for these two equilibration stages as needed. Production runs were then performed in the NPT ensemble with temperature coupling using a Nosé− Hoover thermostat34,35 and pressure coupling using a Parrinello−Rahmam barostat,36,37 both with a relaxation time of 0.5 ps. Unless otherwise stated, the temperature was kept at 298 K and the system pressure at 1 bar. The configuration was constrained by the LINCS algorithm,38 and neighbor lists were updated every 5 time steps in the equilibration runs and every 10 time steps in the production runs. The nonbonded interaction settings were adopted from ref 25.26 Specifically, a twin-range cutoff scheme (0.8 nm for the short-range cutoff and

swell and associate, forming a gel structure that is cross-linked through the unsubstituted regions. All of these mechanisms involve heterogeneity of the substitution as critical to the gelation mechanism, although they disagree on whether it is the more highly or the less highly substituted regions that form the nodes of the gel network. Several studies have probed the structure of the thermoreversible gel. Fairclough et al.14 have shown that methylcellulose gelation is a kinetically arrested phase separation. Both Bodvik et al.15 and Lott et al.16 observed using cryo-TEM that methylcellulose gels form a fibril structure. Bodvik et al. speculated that these fibrils resemble those formed by natural cellulose, with backbones aligned parallel to minimize the exposure of hydrophobic methyl groups. Lott et al. observed that these fibril structures have a very uniform diameter of around 15 nm, which suggests a nonparallel chain configuration since there is no obvious reason that the aggregation of chains into a parallel bundle should stop at a large but definite fibril diameter of 15 nm. Although a considerable amount of experimental effort has been devoted to characterizing both the solvation and the gelation of methylcellulose, computer simulations of cellulosic polymers in water have been rare. To date, cellobiose has been used as a model to study the conformation of cellulose,17−19 while cellulose oligomers of up to nine monomer units in solution have also been simulated.20−22 Simulations up to microsecond scale have been used to study cellulose fibrils and cellulose melts.23−25 However, none of the studies to date have investigated the interaction between multiple cellulosic oligomers in aqueous solution in order to probe the very important yet unsolved gelation mechanism of methylcellulose. Here we present a simulation study probing the solvation behavior of methylcellulose and attempt to understand better B



Article

along with the energy-minimum configuration values derived from X-ray and NMR experiments in the literature. The averaged angles (⟨φ⟩, ⟨ψ⟩) are (−88.9, −129.9) and (−60.7, −126.5) for GROMOS 45A4 and GROMOS 56Acarbo, respectively. For each angle, the average values for the two force fields are comparable to each other and are in good agreement with both the experimental results and previous simulation results for cellobiose.17,19,20,22,24,39 It should be noted that for GROMOS 56Acarbo, some values of φ are positive, indicating a flip from a β linkage to an α linkage. However, further analysis of the area under the distribution curve reveals that the α linkage configurations constitute fewer than 0.5% of the total configurations sampled. We want to explore the effect of the neighboring monomer’s methylation on the distribution of the dihedral angles characterizing the configuration of the hydroxymethyl group (−CH2OH) or the methoxymethyl group (−CH2OCH3) (if O6 is methylated) at the O6 position, namely, ω (∠O5−C5− C6−O6) and ω* (∠C4−C5−C6−O6). In the publications where the detailed parameters of GROMOS 45A4 and GROMOS 56Acarbo were reported, distributions of these two dihedral angles were parametrized from experimentally determined conformer populations of single AGUs26,27 in which the O6 position is unmethylated. There are three ideal configurations of the hydroxymethyl group, where (ω, ω*) have the values (−60°, 60°), (60°, 180°), and (180°, −60°). Because the carbon at the C5 position is tetrahedrally bonded, it is always the case that ω = ω* − 120°. Thus, using the characteristic angle ω alone is sufficient to represent all three configurations. To account for fluctuations of the angles, we use a two-letter code (GG, GT, TG), where the two letters represent the conformations of ω and ω* to describe the three possible configurations. Here “G” represents “gauche” [ω ∈ (−120°, 0°) and ω ∈ (0°, 120°), where the two angles in parentheses are the lower and upper bounds of the range of values of ω] and “T” represents “trans” [ω ∈ (120°, −120°)]. In Table 1 we compare the preferences for the three configurations in each of the two force fields for cellobiose, 6methylcellobiose (6M-cellobiose), and 2,6-dimethylcellobiose (2,6M-cellobiose) at both room temperature (25 °C) and elevated temperature (90°C). The structure of cellobiose has a few subtle yet important differences from the schematic shown in Figure 1. Cellobiose has only two (n = 2) repeat AGUs, M1 on the left and M2 on the right, and all of the substituent groups (R) are hydrogens. The C4 carbon in M1 is connected to a terminal hydroxyl group (−OH) instead of the O4 ether linkage (−O4/O1′−) as shown in Figure 1. Similarly, the C1′ carbon in M2 is also connected to a terminal hydroxyl group (−OH) instead of the O1′ ether linkage (also designated as the O4 ether linkage). Because the two monomers M1 and M2 in cellobiose have terminal −OH groups on opposite sides of the monomer (on the left for M1 and on the right for M2, as depicted in Figure 1), they form different hydrogen bonds and have different conformations as a result of differing chain-end effects on the right side versus the left side of the dimer. Specifically, we expect that the O6 oxygen in M1 can form an intramolecular hydrogen bond with the nearby terminal −OH group attached to the C4 carbon, which is a hydrogen bond that is not available to the O6′ group in M2 because the O6′ group predominantly forms a hydrogen bond with the O2 group. Analogously, for M2 the O2′ oxygen can form a hydrogen bond with the terminal −OH group attached to the C1′ carbon, while its counterpart O2 oxygen in M1 cannot

1.4 nm for the long-range cutoff) was used to handle the nonbonded interaction. The long-range electrostatics beyond 1.4 nm was handled by the reaction-field method with a dielectric constant set at 61 for SPC water. For the simulations with box sizes smaller than 10 nm, a typical simulation run of 20 ns took about 72 h on one or two eight-core Intel Xeon Phi coprocessors with 2 GB of memory per core. The simulations with box sizes of 12, 20, and 22 nm were carried on Kraken and Trestles, which are supercomputers available under the XSEDE computational gateway program (see the Acknowledgment). Each of these simulations of 35 ns used up to 256 cores, and each finished in about 1 week.

■

RESULTS AND DISCUSSION Conformational Preferences of Cellobiose and Methylated Cellobiose. We first study the dimer cellobiose and its seven distinct homodimer methylated derivatives, where there are three distinct derivatives with one methylation site per monomer, three with two, and one with all three sites methylated. These are simulated in water using both the GROMOS 45A4 and GROMOS 56Acarbo force fields to compare and validate the results of these two force fields. GROMOS 45A4 has been validated quite extensively in the past on the basis of simulations of various disaccharides and oligosaccharides, including cellobiose, amylose fragments, and a selection of short oligomers of methylated cellulose.19,20,27,39 The GROMOS 56Acarbo force field, on the other hand, has not been systematically validated for methylated cellulosic dimers or fragments. Even though the performance of the new force field should not deviate dramatically from that of GROMOS 45A4,26 it is still worth validating it before applying it to simulations of longer cellulosic oligomers. For each of the simulated cellobiose derivatives, the distributions of the two characteristic dihedral angles that define the β(1→4) linkage, namely, φ (∠O5−C1−O4′−C4′) and ψ (∠C1−O4′−C4′−C5′), were computed from two 5 ns simulations of one cellobiose molecule using the two GROMOS force fields, and the results are plotted in Figure 2

Figure 2. Distributions of the dihedral angles φ and ψ from simulations of cellobiose in water using both the GROMOS 45A4 and GROMOS 56Acarbo force fields, along with the two pairs of experimentally determined average values ⟨φ⟩ and ⟨ψ⟩ from the literature. φ = ∠O5−C1−O4′−C4′ and ψ = ∠C1−O4′−C4′−C5′, where the primes denote atoms of the monomer M2). C



Article

Table 1. Conformation Distributions of the Hydroxymethyl (−CH2OH) and Methoxymethyl (−CH2OCH3) Groups at the O6 Position in Three Different Cellobiose and Three Different 10-Mer Cellulosic Oligomers at 25 and 90 °C, Along with the Distributions of the Hydroxymethyl Group Determined from Previous Experimental and Simulation Studies of AGUa 25 °C

90 °C

45A4 cellobiose 10-mer cellulose 6M-cellobiose 10-mer 6-MC 2,6M-cellobiose 10-mer 2,6-MC

Ml M2 middle blocks Ml M2 middle blocks Ml M2 middle blocks

GG

GT

52.5 65.4 66.1 68.9 47.9 63.4 62.7 100.0 96.5

47.5 34.6 33.9 30.8 52.0 36.6 37.3 0.0 3.5

56Acarbo TG

monomer

56Acarbo

GT

TG

GG

GT

TG

GG

GT

TG

0.0 19.0 0.0 48.2 0.0 40.4 0.3 30.8 0.1 45.0 0.0 46.9 0.0 16.2 0.0 58.2 0.0 52.9 experiments

68.9 49.4 57.5 67.6 44.7 45.7 83.5 37.4 43.8

12.1 2.3 2.1 1.6 10.4 7.4 0.3 4.4 3.2

31.4 40.5 63.3 51.0 55.3 68.2 50.6 75.5 62.8

68.6 59.5 36.2 48.9 44.3 31.3 48.9 24.2 37.0

0.0 21.4 0.0 44.3 0.5 38.7 0.0 24.2 0.4 39.1 0.4 42.5 0.5 26.2 0.3 45.9 0.2 43.6 simulations

66.2 52.6 57.4 73.4 48.2 50.7 71.8 45.4 51.7

12.4 3.1 3.9 2.4 12.7 6.8 2.0 8.7 4.7

refs 40 and 41 single AGU

45A4

GG

45A427

ref 42

56Acarbo26

GG

GT

TG

GG

GT

TG

GG

GT

TG

GG

GT

TG

53

45

2

41

52

7

57

43

0

37

60

3

All of the conformations are based on the dihedral angle ω (ω = ∠O5−C5−C6−O6) and are tabulated as percentages, where angles near ω ∈ (−120°, 0°), ω ∈ (0°, 120°), and ω ∈ (120°, 240°) [the last of which is equivalent to (−240°, −120°)], correspond to the GG, GT, and TG configurations, respectively. M1 and M2 in cellobiose denote the first (left) and second (right) monomers shown in Figure 1 with n = 2, respectively. The “middle blocks” values for the 10-mers are averages over the two middle blocks (blocks 5 and 6). The experimental conformation distributions of the hydroxymethyl group of single AGU monomers were measured by Nishida and co-workers40,41 and Stenutz et al.,42 and are listed in the first half of the bottom section of the table. The results of Nishida and co-workers were used by Lins and Hünenberger to parametrize the GROMOS 45A4 force field,27 and the results of Stenutz et al. were used by Hansen and Hünenberger to parametrize the 56Acarbo force field;26 the resulting predictions of these two force fields are given in the second half of the bottom section of the table. a

form the corresponding hydrogen bond with O1 since the O2 group mainly forms hydrogen bonds with the O6′ group. We therefore calculate the angle distributions for the O6 and O6′ groups separately for cellobiose. At room temperature, the GROMOS 45A4 force field shows a consistent preference for the GG conformation over the GT conformation for both M1 and M2, except for M2 in 6M-cellobiose. The TG conformation is nearly absent in both monomers. For GROMOS 56Acarbo, on the other hand, in M1 the GT conformation is greatly preferred over the GG conformation in all three cellobiose molecules, while in M2 GT and GG are nearly equally preferred. A small population of the TG conformation (0.3 to 12.1%) for the O6 group is present in all three cellobiose molecules. The TG conformation, although experimentally observed (refer to the bottom section of Table 1), is captured only by the GROMOS 56Acarbo force field because the parameters for the torsional potential of the hydroxymethyl group in this improved force field were reoptimized from those in the GROMOS 45A4 force field. With the refined hydroxymethyl group torsional potential parameters, the simulation results for GROMOS 56Acarbo correctly show a small population of the TG conformation, which allows the formation of the intramolecular O6−O4 hydrogen bond in the left terminal monomer (M1). In GROMOS 56Acarbo, the percentage of TG conformations in M1 declines drastically in the two cellobiose molecules methylated at the O6 position compared with unmethylated cellobiose. In M2, however, the GG and GT conformations are equally probable in cellobiose and 6M-cellobiose and the GG conformation is more probable than the GT conformation in 2,6M-cellobiose, possibly as a result of the interaction between the O6′ group and the O2 group, which is possible only when the hydroxymethyl group at the O6′ position is in the GG or

TG conformation (the TG conformation is rare, leaving the GG conformation as the most probable way that the O6′ and O2 groups can interact). There are differences between the hydroxymethyl conformations of M1 and M2 in all of the dimer molecules simulated using both force fields because of the end-group effect. For the longer oligomers, the end-group effects will be strongest on the terminal monomers, but interior monomers close to the end monomers will also be affected as a result of conformational interactions between neighboring monomers that propagate inward from the end monomers. To minimize the bias from the end groups, we simulated 10-mer oligomers and calculated the hydroxymethyl group conformational preferences of the middle two monomers. Since the conformational preferences for the two neighboring middle monomers are very similar as expected (data not shown), we therefore show in Table 1 only the preferences averaged over these two middle monomers. The majority of the preferences calculated for the middle blocks in the 10-mer oligomers are close to those observed for M2 in the dimer molecules, except for 6-MC with GROMOS 45A4, where the trends for the 10-mer are closer to those observed for M1 in 6M-cellobiose. We have thus shown here that dimers are indeed not good models for determining conformational preferences for the whole oligomer chain. Matthes et al.17 have also suggested that cellobiose is not a good model to study the conformational preferences of cellulose. In the rest of this study, we therefore consider only oligomers that have 10 or more monomer units. Among the three 10-mer cellulosic oligomers studied here, there are differences between the hydroxymethyl conformations obtained using the two GROMOS force fields, which originate from the different experimental data used to parametrize the two force fields. The parameters in GROMOS 56Acarbo were D



Article

obtained by reoptimization of GROMOS 45A4 to make the conformational preferences of the hydroxymethyl group match those obtained from recent experiments on single AGUs (see refs 132−134 in ref 26). In particular, early experimental studies of single AGUs in water, listed as the first set of data in the “experiments” portion of the bottom section of Table 1, under the heading “refs 40 and 41”, showed a preference for the GG conformation over the GT conformation. However, more recent experimental studies show that the GT conformation is preferred over GG, as shown by the second set of data in the “experiments” portion of the bottom section of Table 1 under the heading “ref 42”. The disagreement between these experimental studies is reflected in the different preferences for these hydroxymethyl conformations produced by the two force fields applied to a single AGU, as shown in the “simulations” portion of the bottom section of Table 1 under the headings “45A4” and “56Acarbo”. In our simulations of oligomers, this difference in hydroxymethyl conformational preferences is manifested in the different preferences for the GG or GT conformation in the cellulose oligomers for the two force fields. Nevertheless, at room temperature, GG is preferred over GT in both 6-MC and 2,6-MC methylated cellulose for both force fields. For GROMOS 56Acarbo, when the temperature is increased from 25 to 90 °C, cellulose shows an increase in the probability of the TG conformation, a decrease in the probability of the GG conformation, and an essentially unchanged probability of the GT conformation. These temperature dependences of the conformational preferences are in agreement with a previous simulation study on a cellulose oligomer using the AMBER force field,22 although in that study the GT conformation showed a slightly increased probability at elevated temperature. A similar trend, namely, an increase in the probabilities of the GT and TG conformations and a decrease in the probability of GG at elevated temperature, is shown in Table 1 for methylated cellulose, except for the TG conformation in 6-MC. Persistence Length and Chain Length Dependence of Rg. The persistence length has been measured experimentally for methylcellulose but has not been estimated from previous simulations because of limitations on the oligomer chain lengths that could be simulated. Figure 3a shows our calculations of the radius of gyration (Rg) as a function of chain length and the fit of the persistence length from the Kratky−Porod wormlike chain model for MC oligomer chains, both obtained using the GROMOS 45A4 force field (results obtained using the 56Acarbo force field are shown in Figure SI1 in the Supporting Information). Two sets of MC oligomers are simulated, namely, 2,6-MC homo-oligomers and random MC co-oligomers. The 2,6-MC homo-oligomers with 10, 20, 35, and 40 monomer units are simulated using GROMOS 45A4, while homo-oligomers with 8, 15, 25, 35, and 40 repeat units are simulated using 56Acarbo. The two random co-oligomers with 18 and 28 repeat units are simulated using both force fields. The 2,6-MC homo-oligomer is chosen because of the high abundance of this monomer in the METHOCEL A chemistry. The oligomers are initially in the all-trans configuration in all of the simulations and are allowed to relax for 10 ns (for oligomers with chain lengths less than 20) or 20 ns (for the rest of the oligomers) before the data used to compute Rg are gathered. The length of each cubic simulation box is set equal to the contour length of the oligomer, which is estimated to be L = 0.54n nm, where n is the number of monomers in the chain and the factor of 0.54 was reported by

Figure 3. (a) Simulated radii of gyration (Rg) of MC chains obtained using the GROMOS 45A4 force field and fits to the Kratky−Porod model using various persistence lengths lp. (The fit of simulated Rg values obtained using GROMOS 56Acarbo is shown in Figure SI1 in the Supporting Information.) (b) Rg values of 2,6-MC and random MC, where results for all of the molecules simulated using both force fields are plotted against molecular weight along with experimental data. (c) Lengths of cubic simulation boxes used for the single-chain simulations.

Patel et al.43 The radii of gyration (Rg) are computed using the g_gyrate function in GROMACS and averaged over two 5 ns periods. If the averaged Rg values obtained from the two 5 ns intervals differ by more than 10%, the oligomer is simulated for another 5 ns to ensure the convergence of the averaged Rg values. The simulated Rg values for all of the oligomers obtained using GROMOS 45A4 are plotted against those for a rodlike oligomer in Figure 3a, where Rg is estimated to be L/ √12. The data from the two force fields are very similar and start to deviate from rodlike behavior beyond 20 monomers length, which roughly corresponds to the persistence length of MC chains. E



Article

Table 2. Simulated Average Dihedral Angles (in deg) with Standard Errors for Cellulose and Methylated Cellulose at 25 and 90 °C Obtained Using the GROMOS 45A4 and 56Acarbo Force Fieldsa 25 °C

90 °C

45A4 cellulose 2-MC 3-MC 6-MC 2,3-MC 2,6-MC 3,6-MC 2,3,6-MC a

56Acarbo

45A4

56Acarbo

⟨φ⟩

⟨ψ⟩

⟨φ⟩

⟨ψ⟩

⟨φ⟩

⟨ψ⟩

⟨φ⟩

⟨ψ⟩

−71.8(3) −72.8(2) −76.1(4) −73.3(2) −77.2(2) −73.0(2) −76.4(2) −105.4(2)

−126.4(2) −125.1(1) −115.7(2) −129.1(1) −104.4(3) −128.6(1) −89.0(2) −126.2(3)

−64.6(3) −66.2(2) −57.5(3) −66.4(2) −60.4(3) −66.7(2) −49.7(3) −64.4(2)

−128.9(2) −127.4(1) −125.9(2) −134.4(2) −125.6(2) −128.5(1) −122.5(1) −105.2(2)

−71.5(2) −71.5(3) −81.8(3) −70.4(2) −81.4(2) −70.9(2) −76.6(2) −95.7(3)

−126.1(1) −119.7(2) −114.8(4) −120.7(2) −107.8(5) −127.5(1) −102.9(3) −112.0(4)

−70.9(4) −67.8(2) −53.0(3) −59.4(3) −65.9(2) −65.7(2) −56.8(4) −71.3(2)

−126.5(5) −128.5(2) −126.3(2) −128.1(2) −129.2(2) −128.5(2) −123.6(2) −99.7(4)

In this and the following, the error in the last digit is indicated in parentheses, e.g., −71.8(3) indicates −71.8 ± 0.3.

We further fit the simulated Rg values with the Kratky−Porod model44 for semiflexible chains, for which: Rg2 =

2lp 4 2lp3 1 lpL − lp2 + 2 (e−L / l p − 1) + 3 L L

estimated persistence length. However, we believe that both force fields are able to reproduce the persistence length of methylcellulose chains quite well. Although there are no experimental values for Rg at molecular weights as low as those of our simulations, in Figure 3b we plot our results along with the experimentally determined Rg values for various MC polymers with molecular weights in the range 20−1200 kDa from both Keary5 and Patel et al.43 The Rg values reported by Patel et al. fall in the range between the upper and lower Rg trend bounds predicted using the Kratky−Porod model based on the obtained persistence length estimation from the simulations by both force fields (7.6−15.5 nm). The Rg values reported by Keary, however, are slightly below the lower predicted Rg value trend line bound (lp = 7.6 nm). A fit of Keary’s data using the Kratky−Porod model yields a persistence length of around 6 nm (Figure 3b inset), which is lower than the persistence length reported in the literature and estimated from our simulation. It should be noted that Keary pointed out that the Flory−Fox equation used to estimate Rg may underestimate the Rg values of cellulosic chains because the equation is applicable to a linear flexible polymer chain while cellulosic chains are relatively rigid. Three regimes can be observed in Figure 3b: (1) the rodlike regime at very low molecular weight, which is observed only in the simulation data; (2) the flexible regime, where the contour length is much higher than the persistence length, which is observed only in the experimental data; and (3) the transition regime between the two, which lies roughly in the range of molecular weights in the gap between the simulation data and the experimental data. Fitting the experimental Rg values reported by Keary to a power-law expression yields an exponent of 0.57, which corresponds to the Rg scaling for flexible polymers in a good solvent, for which the exponent should be around 0.59. The experimental Rg values at the lower end of the molecular weight range seem to match up well with the Rg values obtained for 35mers and 40-mers from the simulations, indicating that a chain length of around 40 monomers is at the crossover between the stiff- and flexible-polymer regimes. It should be noted that for a semiflexible polymer such as methylcellulose, the Kuhn length is twice the persistence length. The Kuhn length estimated from our simulation is around 20 nm, which corresponds to a 40-mer chain. Solvation Behavior of a Single Oligomer Chain. We next explore the solvation behavior of single MC chains. Although we are able to simulate a single 40-mer model oligomer, doing so requires significant computational effort and

(1)

The persistence length (lp) can be estimated by fitting this model to the simulated Rg values of oligomers with different chain contour lengths (L). The fit was conducted using the curve-fitting toolbox in MATLAB, version R2012a (The MathWorks). The persistence length estimated for model oligomers simulated using the GROMOS 45A4 force field is 11.9 nm, with a 95% confidence interval (CI) of 8.4−15.5 nm. This estimation is in good agreement with the experimental persistence lengths obtained for six different methylcellulose samples by Patel et al., ranging between 12 and 17 nm with error bars of ±2 nm for each sample.43 The estimated persistence length for model oligomers simulated using GROMOS 56Acarbo is 9.4 nm, with a 95% CI from 7.6 to 11.3 nm. This value is also in reasonable agreement with the reported experimental persistence length of MC. The estimated persistence lengths obtained from the 56Acarbo force field are somewhat smaller than those from the 45A4 force field, although the confidence intervals of the two overlap. We suspect that this difference, if it is real despite the overlapping confidence intervals, could be due to the following. First, the dominance of the GG conformation of the hydroxymethyl group in the GROMOS 45A4 force field allows a strong hydrophobic interaction between O6 and O2′, or an intramolecular O6−O2′ hydrogen bond if either of the two groups is not substituted (an O6−O2′ hydrogen bond is equivalent to an O6′−O2 hydrogen bond), which helps to maintain the chain stiffness. This chain-stiffness stabilization is weakened in the GROMOS 56Acarbo force field because half of the time the hydroxymethyl group is in the GT conformation, which allows the O6 group to interact only with the O3′ group. The distance between the O6 and O3′ groups is greater than that between the O6 and O2′ groups; thus, neither a hydrophobic interaction nor a hydrogen bond is likely to be common between O6 and O3′ groups (O6 and O6′ are chemically identical, and the prime indicates that the group is on the right monomer of two consecutive monomer units on the cellulosic chain). Although it is possible that the two force fields give different persistence lengths for the above reason, it should be noted that the estimated persistence length is fairly sensitive to the simulated Rg value, and a small difference in the Rg value (∼0.4 nm) could lead to a 4 nm difference in the F



Article

thus is not very cost-effective. We therefore instead simulate 10mer cellulose and methylated cellulose chains using the two GROMOS force fields at both room temperature and elevated temperature (90 °C) to study the effect on the solvation behavior of different substitution patterns and temperatures. Each system is equilibrated for 10 ns before the data are collected. Table 2 summarizes the simulated average dihedral angles of these oligomers under all simulation conditions. Each value is the average over both the linkage between the fourth and fifth monomers and that between the fifth and sixth monomers and is also averaged over the simulation time after equilibration of the chain configuration. The results obtained from both force fields show the same trends. At room temperature, the average dihedral angles ⟨φ⟩ and ⟨ψ⟩ of the oligomers that do not have the O3 position substituted (i.e., cellulose and 2-, 6-, and 2,6MC) are very similar. (Hereafter, we leave off the brackets “⟨ ⟩” denoting that the angles are averages.) However, if the O3 position is substituted, some dihedral angles, namely, ψ for the methylated oligomers simulated using GROMOS 45A4 and φ for those simulated using GROMOS 56Acarbo, deviate from those of natural cellulose, with the exception of 2,3,6-MC, where φ for oligomers simulated using GROMOS 45A4 and ψ of those simulated using GROMOS 56Acarbo deviate from the corresponding dihedral angles of natural cellulose. These deviations occur because methylcellulose oligomers with the O3 position substituted (i.e., 3-, 2,3-, and 2,3,6-MC) lack the O3′−O5 hydrogen bond (as shown in Figure 1, the O3′−O5 hydrogen bond is equivalent to the O3−O5′ hydrogen bond), which is the hydrogen bond that helps to maintain the planar geometry of the oligomer. At elevated temperature, the dihedral angles do not change dramatically from those at room temperature, which shows that heating to 90° does not cause additional bending or distortion of the oligomers, at least for chains of 10 monomers. Yu et al.20 have reported that the dihedral angles of selected single-chain 9-mer cellulosic oligomers, namely, cellulose and 6-, 2,3-, and 2,3,6-MC oligomers, simulated using the GROMOS 45A4 force field in water are (φ, ψ) = (−94°, −134°), (−99°, −137°), (−101°, −139°), and (−104°, −143°), respectively. The corresponding results in our study agree with these values to within around 20°, with the exception of ψ in 2,3-MC, where a 35° difference is found. These differences possibly originate from the longer equilibration time used in our study (10 ns) compared with that used by Yu et al. (500 ps). We compute the oligomer−water radial distribution functions (RDFs) of cellulose and methylated cellulose molecules at both 25 and 90 °C using the GROMACS function g_rdf, and the results at 25 °C for both GROMOS force fields are plotted in Figure 4. The oligomer−water RDF is the density of atoms on water molecules at a distance r from an atom on the oligomer, averaged over all atoms of the oligomer and over time. The RDF is cut off at radius rmax (half of the box side length, 3 nm) as the value of g(r) approaches 1. There are two key features in these oligomer−water RDF curves: the height of the RDF curve at 1 nm (main graph) and the height of the nearest-neighbor peak between 0.15 and 0.20 nm (which is cut off in the main figure but enlarged in the insets of Figure 4). The heights of the RDFs in the main graphs correlate well with the number of oligomer−water hydrogen bonds (see the next paragraph for the definition of hydrogen bonds) formed by cellulose and MCs with different degrees of substitution, which are shown in Table 3. It is not surprising that cellulose has the

Figure 4. Oligomer−water radial distribution functions (RDFs) of 10mer cellulose and methylated celluloses at 25 °C obtained using (a) GROMOS 45A4 and (b) GROMOS 56Acarbo. The insets are the zoomed-in views near the first peak of the RDF curves (between r = 0.1 and 0.25 nm). The ordering of the legend, from top to bottom, matches the heights of the curves, in both the main figures and the insets.

highest RDF peak, as all of the hydroxyl groups remain unsubstituted and are therefore available to form hydrogen bonds (since only one short chain is present, the majority of these hydrogen bonds are formed with water molecules rather than with other cellulose molecules, as occurs in cellulose crystals). The monosubstituted methylcellulose molecules have lower RDF peak heights, followed by the disubstituted ones and then the trisubstituted one. The nearest-neighbor peaks also correlate with the number of hydrogen bonds formed between the oligomer and the solvent, with clear differences among methylcellulose oligomers with the same DS, namely, the three monosubstituted MCs and the three disubstituted MCs (Table 3). Of the monosubstituted MCs, 3-MC has the highest peak, followed by 2-MC and finally 6-MC. This suggests that the three hydroxyl groups have different hydrogen-bond-forming propensities. Thus, methylating the O3 oxygen has less effect on the number of oligomer−water hydrogen bonds formed than does methylating the other oxygens. The same conclusion holds for the three dimethylated cellulose molecules, where 2,6-MC has a lower nearestG



Article

Table 3. Average Numbers of Intermolecular (Oligomer−Water, O−W) and Intramolecular (Oligomer−Oligomer, O−O) Hydrogen Bonds with Standard Errors for a Single 10-Mer in Water at 25°C and 90°C Using the Two GROMOS Force Fields 45A4

56Acarbo

25 °C cellulose 2-MC 3-MC 6-MC 2,3-MC 2,6-MC 3,6-MC 2,3,6-MC

90 °C

25 °C

90 °C

O−W

O−O

O−W

O−O

O−W

O−O

O−W

O−O

61.87(6) 47.91(5) 51.53(5) 46.84(5) 38.65(5) 32.93(4) 40.33(5) 28.74(4)

5.30(2) 5.65(2) 1.70(2) 7.50(2) 0.85(1) 7.98(1) 0.27(1) 0.00(0)

53.75(6) 42.16(6) 44.11(6) 40.30(6) 32.78(5) 28.72(5) 34.89(5) 23.80(5)

5.86(2) 5.76(2) 1.98(2) 7.30(2) 1.29(1) 7.39(2) 0.54(1) 0.00(0)

60.96(6) 46.93(5) 50.17(6) 46.68(5) 37.96(5) 32.79(5) 41.53(5) 29.24(4)

5.30(2) 6.07(2) 2.97(2) 6.97(2) 2.73(2) 7.34(2) 0.37(1) 0.00(0)

52.38(6) 40.76(6) 43.33(6) 40.38(6) 31.66(6) 28.07(5) 35.27(5) 23.09(5)

5.55(2) 6.20(2) 3.06(2) 6.26(2) 2.75(2) 7.12(2) 0.65(1) 0.00(0)

Figure 5. Snapshots of various methylcellulose oligomers in water at 25 °C after 35 ns of simulation time.

intermolecular hydrogen bonds formed between oligomers and water and the number of intramolecular hydrogen bonds formed within the single oligomer using the GROMACS function g_hbond with the definition mentioned above, and the results are summarized in Table 3. The numbers of intermolecular and intramolecular hydrogen bonds obtained using the the two force fields agree very well. We have already shown that the numbers of intermolecular oligomer−water hydrogen bonds correlate well with the oligomer−water RDF peak heights. In addition, intramolecular hydrogen bonding has been shown to stabilize intramolecular cellulosic structures.19−22 A detailed analysis of the hydrogen bonds formed

neighbor peak than the other two, again suggesting that methylation of the O3 oxygen has less effect on the number of oligomer−water hydrogen bonds formed than does methylation of the O2 and O6 oxygens. This conclusion is in agreement with a recent simulation study on cellobiose.19 The oligomer−water RDFs at 90 °C show the same trend with negligible shifts of the peak positions and thus are not shown here. In this study, a hydrogen bond is defined to exist when (1) the distance between the donor and acceptor oxygen atoms is less than 3.5 Å and (2) the hydrogen atom−donor−acceptor angle is smaller than 30°.45,46 We compute the number of H



Article

Figure 6. Same as Figure 5 except at 90 °C.

Snapshots of the final configurations at 35 ns for cellulose and methylcellulose oligomers in water at 25 and 90 °C are shown in Figures 5 and 6, respectively. At 25 °C, cellulose and monosubstituted methylcellulose show loosely aggregated structures, while di- and trisubstituted methylcelluloses show more densely aggregated clusters. At 90 °C, all of the oligomers except cellulose, 3-MC, and 6-MC form what appear to be precipitates. The oligomers depicted in Figures 5 and 6 are homo-oligomers, which may have different solvation behavior than the experimental random hetero-oligomers. Experimentally, at room temperature methylcellulose is only water-soluble for DS between roughly 1.4 and 2.0. At 90 °C, which is well above their experimentally determined gelation temperatures, methylcelluloses typically exhibit precipitate-like structures.14,15 In what follows, we use oligomer−oligomer RDFs and numbers of hydrogen bonds to further characterize these multiple oligomer systems. We compute oligomer−oligomer RDFs by first finding all atoms on any oligomers at a distance r from an atom on a given oligomer and excluding any atoms that connect to the starting atom through a shortest sequence of 15 or fewer bonds. We then compute the density of atoms on the basis of this criterion and average over all of the oligomer atoms and then average over time. The resulting RDFs are insensitive to the bond exclusion number when it is larger than 15, indicating that the contributions to the RDF from nearby chemically bonded monomers are fully excluded using this bond exclusion number. The RDF is cut off at radius rmax (half of the box side length, 6 nm) as g(r) approaches 1. The RDFs are computed using the

at each oxygen position reveals that the O3 position mostly forms the O3−O5′ intramolecular hydrogen bond, which is one of the predominant hydrogen bonds that helps to maintain the planar geometry of the oligomer (the O3−O5′ hydrogen bond depicted in Figure 1 is equivalent to the O3′−O5 hydrogen bond). We also find that at elevated temperature the simulation box expands 6% by volume (data not shown) at constant pressure as water comes close to its boiling point. This reduces the number of water molecules that are close to the oligomer, and as a result, the oligomer forms fewer oligomer− water hydrogen bonds. The number of intramolecular hydrogen bonds is not altered when the simulation temperature is raised (except for the case of cellulose). This is expected because the characteristic dihedral angles φ and ψ are not much affected by the temperature. Simulations Containing Multiple Homo-oligomers. Here we present results for the solvation behavior of multiple oligomers in water. Because the two GROMOS force fields yield similar molecular conformations and numbers of hydrogen bonds, we limit further simulations to the GROMOS 56Acarbo force field with short 10-mers for computational efficiency. Each simulation is carried out in a 12 nm cubic simulation box with 27 oligomers initially randomly distributed in the box, bringing the oligomer concentration to 4.5 wt %. Simulations are carried out at two temperatures, 25 and 90 °C, to explore the temperature dependence. To our knowledge, this is the first attempt to use atomistic simulations to study the solvation behavior of multiple (>20) methylcellulose chains and to explore temperature effects. I



Article

aggregated structure at high temperature similar to that of disubstituted methylcellulose. We believe that this is related to the formation of intermolecular hydrogen bonds in 2-MC at high temperature, which is further discussed in the following. We compute the number of hydrogen bonds and categorize them as oligomer−water, intramolecular oligomer−oligomer, and intermolecular oligomer−oligomer hydrogen bonds. The total numbers of oligomer−water and oligomer−oligomer hydrogen bonds are divided by the number of chains to normalize for system size in Table 4. The number of oligomer−

g_rdf function implemented in GROMACS, and the RDF curves themselves are shown in Figure SI2 in the Supporting Information. We find that the peak of the RDF (which occurs at roughly r = 1.2 nm) levels off as a function of time after 30 ns, suggesting that a relatively stable oligomer aggregate has formed (see Figure SI3 in the Supporting Information). As a quantity that differentiates visually similar structures, especially among MCs with the same DS at the same temperature, in Figure 7 we plot the RDF peak heights of all eight homo-

Table 4. Total Numbers of Oligomer−Water (O−W) and Oligomer−Oligomer (O−O) Hydrogen Bonds, Normalized by the Number of Chains, in Various Oligomers at 4.5 wt %, with Standard Errors 25 °C cellulose 2-MC 3-MC 6-MC 2,3-MC 2,6-MC 3,6-MC 2,3,6-MC

90 °C

O−W

O−O

O−W

O−O

53.6(9) 41.6(7) 43.6(7) 41.1(7) 27.8(6) 23.1(5) 35.1(7) 17.4(5)

7.8(4) 7.3(3) 5.1(3) 8.1(3) 4.0(2) 8.2(2) 1.5(2) N/A

48.4(9) 31.2(9) 35.2(9) 33.2(9) 19.5(6) 16.3(5) 24.8(6) 10.5(5)

6.9(4) 8.4(4) 5.7(4) 8.1(3) 4.4(3) 8.1(2) 2.1(2) 0.2(1)

water hydrogen bonds follows the same trend as in the corresponding single-oligomer cases, and these are the most abundant hydrogen bonds formed. Table 5 shows a detailed breakdown of the types and numbers of oligomer−oligomer hydrogen bonds per oligomer chain. The majority of the oligomer−oligomer hydrogen bonds are intramolecular ones, ranging from 42 to 97%. However, there are also many intermolecular hydrogen bonds. At the elevated temperature (90 °C), there is an increase in the percentage of intermolecular hydrogen bonds for most of the cases (except for cellulose and 3,6-MC), which we believe is partially responsible for the kinetically arrested phase separation (gelation) of the methylated polymer at high temperature. We further distinguish the intermolecular hydrogen bonds on the basis of the three different substitution positions. We find that the O6 position, if not methylated, forms many intermolecular hydrogen bonds with other unmethylated substitution positions. In particular, the O6−O6 intermolecular hydrogen bond is the most abundant intermolecular hydrogen bond in most systems if its formation is possible (e.g., in cellulose, 2-MC, and 2,3-MC; also see Table SI1 in the Supporting Information for a complete summary of intermolecular hydrogen bond percentages). O2−O2 and O2−O6 are the next most probable, presumably because the O2 position has a hydrogen-bondforming tendency that is second to that of the O6 position. It should be noted that one exception to this observation occurs in the 3-MC complex, where the most probable intermolecular hydrogen bonds are O2−O2 and O2−O6 at room and elevated temperature, respectively. This is surprising because the O6− O6 intermolecular hydrogen bond is also possible in this complex and yet occurs less frequently than the two possible intermolecular hydrogen bonds formed at the O2 position (i.e., O2−O2 and O2−O6). We speculate that the presence of the methylated O3′ position on the neighboring monomer somehow interferes with intermolecular hydrogen-bond formation at the O6 position.

Figure 7. Heights of the peaks of the oligomer−oligomer RDFs (at r ≈ 1.2 nm) for the various oligomers in water plotted against the number of oligomer−water hydrogen bonds in the corresponding systems. The blue and red data points are for RDF peaks at 25 and 90 °C, respectively. The legend is ordered from the highest RDF peak height to the lowest at 25 °C.

oligomer systems at two different temperatures against the numbers of oligomer−water hydrogen bonds formed at 25 °C (although the number of oligomer−water hydrogen bonds depends on temperature, we plot against the number of oligomer−water hydrogen bonds at 25 °C on the x axis to line up vertically the results for each oligomer type at the two temperatures). On the basis of this metric, at 25 °C 2-MC has the “loosest” structure, with the lowest peak height, followed by cellulose, 6-MC, 3-MC, and then the di- and trisubstituted methylcelluloses. This is expected since disubstituted and trisubstituted methylcellulose molecules have more hydrophobic methyl substituent groups than the monosubstituted methylcellulose molecules and thus have a greater tendency to cluster. At 90 °C, each RDF peak is higher than that at 25 °C except for cellulose, which has a lower RDF peak than at low temperature. Cellulose dissolves better at elevated temperature because of the decreased number of intermolecular oligomer− oligomer hydrogen bonds formed. Although cellulose is expected to crystallize, neither the size of the simulation box nor the duration of the simulation permit this. Various methylcellulose systems show very tight structures at 90 °C, with 2,3,6-MC and 2,3-MC having the highest peak heights, followed by 2,6-MC, 3,6-MC, and 2-MC. The increase in the RDF peak heights with temperature for all of the methylcellulose solutions clearly demonstrates that these simulations are able to reflect the temperature-dependent solvation behavior and hints at aggregation or even precipitation at high temperature (>75 °C). It is surprising that 2-MC forms an J



Article

Table 5. Percentages of Various Kinds of Oligomer−Oligomer Hydrogen Bonds; O2, O3, and O6 Represent the Intermolecular Hydrogen Bonds Formed at the Respective Substitution Positionsa 25 °C cellulose 2-MC 3-MC 6-MC 2,3-MC 2,6-MC 3,6-MC 2,3,6-MC a

most probable intermolecular hydrogen bond

90 °C

intra

inter

O2

O3

O6

intra

inter

O2

O3

O6

25 °C

90 °C

68 85 71 91 85 97 42 >99

32 15 29 9 15 3 58

40 2 52 56 5

16 31 1 32 7 6 1

44 67 46 12 88 94 11

80 76 63 88 76 96 50 >99

20 24 37 12 24 4 50

35 4 50 63 3

17 33 2 29 2 15 3

48 63 48 8 95 85 12

O6−O6 O6−O6 O2−O2 O2−O2 O6−O6 N/A O2−O2 N/A

O6−O6 O6−O6 O2−O6 O2−O2 O6−O6 N/A O2−O2 N/A

88

84

Empty entries indicate that no intermolecular hydrogen bonds formed for the particular configurations.

It is also noteworthy that the number of oligomer−water hydrogen bonds does not inversely correlate with the height of the oligomer−oligomer RDF peak. Specifically, in both singlechain and multiple-chain simulations, the 3,6-MC oligomer forms the highest number of oligomer−water hydrogen bonds among the disubstituted MCs, followed by the 2,3-MC oligomer and finally the 2,6-MC oligomer, as shown in Table 4. However, for multiple-oligomer simulations, 2,3-MC has the highest oligomer−oligomer RDF peak at both temperatures, followed by 2,6-MC and finally 3,6-MC, as shown in Figure 7. This observation again indicates that intermolecular hydrogen bonds are partially responsible for the different solvation behaviors of methylcellulose molecules, despite similar degrees of substitution (e.g., DS = 2). Simulation of Multiple Random Oligomers Mimicking METHOCEL A Chemistry. We have shown that MD simulations using the GROMOS 56Acarbo force field are able to discriminate the effects of different methyl-group substitution patterns on the peaks of RDF curves and on the numbers of hydrogen bonds for homo-oligomers of MC. We now apply this force field to simulate model hetero-oligomers whose monomers occur with probability matching that of the METHOCEL A chemistry. To do so, we set up four model systems with oligomer lengths of 10, 20, and 40 monomers. To increase the randomness of our model oligomers, we first generate one 90-monomer random sequence and break it into nine 10-mers, each with a different sequence. Similarly, we generate two 400-monomer sequences and break one of them into 20 random 20-mers and the other into 10 random 40mers. Each initial sequence of 90 or 400 monomers contains monomer probabilities based on the monomer probabilities of METHCOEL A chemistry, defined in Table 6, with a DS of

around 1.9. We choose 20 nm cubic boxes for the 20-mer systems and 22 nm cubic boxes for the 40-mer system so that the side of each box is larger than both the length of the oligomer chain and the diameter of a methylcellulose fibril gel reported in the literature (15 nm).16 (While we initially had some hopes of seeing evidence of such fibrils in our simulations, we were not able to detect them, probably because of inadequate length and time scales.) The oligomer concentrations range from 3 to 6.1 wt %. Specifically, the simulations labeled “10-mer 4 wt %”, “20-mer 3 wt %”, “20-mer 6 wt %”, and “40-mer 4 wt %” contain the oligomers at concentrations of 4.5, 3.1, 6.1, and 4.4 wt %, respectively. Although these concentrations are slightly higher than that in the solutions used for many experiments (2 wt %), these higher concentrations are chosen to increase the number of oligomers in boxes that are small enough to be simulated with available computational power. Nevertheless, the experimentally determined gelation temperature is known to be between 40 and 45 °C for methylcellulose solutions in the range of concentrations simulated here. We therefore perform simulations at 25 and 50 °C, which are below and above the experimental gelation temperature, respectively. Snapshots of the four systems at 25 and 50 °C are shown in Figures 8 and 9, respectively. We observe that the “10-mer 4 wt %” system, which has a DS of 1.87, behaves similarly to the 10mer multiple homo-oligomer disubstituted methylcellulose systems (DS = 2) that we simulated earlier, with chains forming a globular aggregate at both the low and high temperatures. The “20-mer 3 wt %” system, however, forms a ramified aggregate at low temperature and a more rodlike structure that spans the box at high temperature. The “20-mer 6 wt %” and “40-mer 4 wt %” systems form three-dimensional “gel-like” networks at both low and high temperatures. At the elevated temperature in Figure 9, the snapshots for these two systems do not appear significantly different from those at low temperature, suggesting that these two systems might be trapped in a local energy minimum. In Figure 10 we show zoomed-in views of the “40-mer 4 wt %” and “20-mer 6 wt %” structures at 50 °C, which reveal that the oligomers condense into bundles with parallel chain alignment, seemingly to minimize the contact between the hydrophobic methyl groups and water. Although this parallel configuration of polymer associating chains agrees with the proposed structure-formation mechanism by Bodvik et al.,15 who observed a methylcellulose fibril structure in cryo-TEM images, we suspect that this is not the actual configuration the flexible methylcellulose chains adopt in the cross-linked gel structure for the following reasons.

Table 6. Average Mole Fraction of Each Methylated Cellulose Monomer in METHOCEL A Chemistry average mole fraction cellulose 2-MC 3-MC 6-MC 2,3-MC 2,6-MC 3,6-MC 2,3,6-MC total DS(Me)

0.05 0.13 0.02 0.10 0.10 0.26 0.05 0.29 1.94 K



Article

Figure 8. Snapshots after 35 ns for simulations of METHOCEL A in water at 25 °C.

Figure 10. Zoomed-in views of the METHOCEL A “40-mer 4 wt %” (top) and “20-mer 6 wt %” (bottom) systems at 50 °C, showing that model METHOCEL A oligomers condense parallel to each other. Each chain is colored differently. Figure 9. Same as Figure 8, except at 50 °C.

Thus, in our simulations the minimization of methyl-group hydrophobic interactions and the formation of intermolecular hydrogen bonds occurs most readily through parallel configurations of different short chains rather than the formation of a helical structure by a single or a few much longer chains. We suggest that although the parallel chain configurations observed in our short-chain simulations may not be representative of the final gel structures that much longer chains would form, they might be precursors of the final selfassembly of methylcellulose chains at larger length scales and longer time scales, which are only accessible through coarsegrained simulations. Srinivas and Klein47 have reported such self-reorganization from clusters of parallel chains to helically twisted columns in their coarse-grained simulations of amphiphilic molecules, and something similar might occur

First, the parallel configuration observed in our simulations does not explain the observation of fibrils with a uniform diameter of 15 nm in methylcellulose, as reported by Lott et al.,16 because there is no apparent reason that a bundle formed by chains condensing parallel to each other would thicken to a specific diameter of 15 nm rather than continuing to thicken indefinitely. Second, the flexible methylcellulose chains in METHOCEL A products have chain lengths ranging from hundreds to thousands of monomer units, which allows the formation of helices with radii in the nanometer range (i.e., 15 nm) as a result of self-attracting forces. The model oligomers simulated here (10-mer, 20-mer, and 40-mer) are too short, given their persistence lengths, to bend enough to form helices. L



Article

4 wt %” and “40-mer 4 wt %” simulations are in boxes of different sizes, namely, 12 and 22 nm, they have very similar RDF peak positions (r ≈ 1.2 nm) and similar peak heights at both low and elevated temperatures. This suggests that for these relatively short oligomers (10 to 40 monomers) the peak height is more sensitive to concentration than to oligomer size and that we can obtain representative results using cheap simulations of “10-mer” oligomers in a small simulation box. We next compute the numbers of hydrogen bonds formed between oligomers and water and also between oligomers and oligomers (Table 7). We again normalize the total hydrogen bond counts by the numbers of 10-mer lengths in each system to allow not only comparison of the results among the four model random co-oligomers but also comparison of the results for the random co-oligomers with those for simulations of multiple homo-oligomer chains reported in the previous section. All four co-oligomers (DS ≈ 1.9) show very similar numbers of oligomer−water and oligomer−oligomer hydrogen bonds. The number of hydrogen bonds at 25 °C, around 29, also agrees nicely with the number in Table 4 at 25 °C, averaged over the three homodimethylated methylcellulose chains (DS = 2). Around 90% of the oligomer−oligomer hydrogen bonds in the co-oligomers are intramolecular (Table 7), and the rest are intermolecular (i.e., between different chains). For all four systems, as the temperature increases, the number of oligomer−water hydrogen bonds decreases and the number of oligomer−oligomer intermolecular hydrogen bonds increases slightly. Thus, similar to the RDFs, the shortest oligomer, namely, a 10-mer, yields the same normalized number of hydrogen bonds and the same temperature dependence of the hydrogen bonding as in the simulations in a much bigger box with longer model oligomers, namely, the “20-mer 3 wt %”, “20-mer 6 wt %”, and “40-mer 4 wt %” systems. Nevertheless, the larger aggregates in the “20-mer 3 wt %”, “20-mer 6 wt %”, and “40-mer 4 wt %” systems are more ramified than those in the “10-mer 4 wt %” solutions, which may be crucial to the gel-formation mechanism, although the simulation length and time scales are nowhere close to those for the experimental gels. Our results indicate that while we are not able to simulate the large-scale structures formed by these gels, the numbers of hydrogen bonds, RDF peaks, and other metrics of local structure on the scale of a few nanometers are insensitive to the length of the oligomers and therefore may be correctly predicted by the simulations. To determine whether methylcellulose gelation is induced by the more hydrophobic monomers or by the more hydrophilic ones in heterogeneous methylcellulose, such as that used commercially, we compute contact maps of our model systems. We define two monomers to be “in contact” if any one atom on

with methylcellulose chains. An interesting experimental test of the mechanism of fibril formation therefore might be to study the gelation of much shorter methylcellulose chains than have been studied heretofore, namely, chains of only 10−40 monomers similar to those studied in our simulations. If short methylcellulose chains such as these do not form 15 nm fibrils but instead form much less regular gel structures, this might support the suggestion here that a minimum chain length is required for the formation of the regular fibrils seen experimentally. We compute the oligomer−oligomer RDFs of all four systems at both room and elevated temperatures, and the heights of the RDF peaks that occur at roughly r = 1.2 nm are plotted against oligomer weight percentage in Figure 11.

Figure 11. Peak heights of oligomer−oligomer RDFs (at r ≈ 1.2 nm) in four METHOCEL A model oligomer solutions at room temperature (25 °C), shown in blue, and elevated temperature (50 °C), shown in red. The oligomer−oligomer RDF is defined earlier in the text. The RDF curves are shown in Figure SI4 in the Supporting Information.

Compared with those at room temperature, the peaks of the RDFs at elevated temperature for all four systems are consistently higher, corresponding to closer-packed aggregates, which is consistent with our earlier observations for homooligomers. In addition, it is clear that the heights of the RDF peaks decrease as the oligomer weight percentage increases. The “20-mer 6 wt %” system, for example, has the lowest RDF peak and shows barely any peak height change in response to the temperature change. Interestingly, even though the “10-mer

Table 7. Numbers of Hydrogen Bonds Per 10 Monomers in the Chain, with Standard Errors, and Percentages of Different Oligomer−OligomerHydrogen-Bond Types in the METHOCEL A Model Systemsa numbers of hydrogen bonds 25 °C 10-mer 20-mer 20-mer 40-mer

4 3 6 4

wt wt wt wt

% % % %

percentages of O−O hydrogen-bond types 50 °C

25 °C

50 °C

O−W

O−O

O−W

O−O

intra

inter

intra

inter

30.7(7) 28.8(5) 28.2(4) 29.7(6)

5.1(2) 5.0(1) 5.2(1) 5.2(1)

26.5(5) 25.3(7) 24.8(7) 25.2(3)

5.1(2) 5.1(1) 5.2(1) 5.4(1)

91 92 92 93

9 8 8 7

89 91 90 92

11 9 10 8

a

The normalization of the hydrogen bond counts per numbers of 10-mer subchains means that numbers of hydrogen bonds in 20-mer and 40-mer chains are obtained from the entries in the table by multiplying by 2 and 4, respectively. M



Article

Figure 12. Contact maps averaged over four model METHOCEL A systems (“10-mer 4 wt %”, “20-mer 3 wt %”, “20-mer 6 wt %”, and “40-mer 4 wt %”) at (a) 25 °Cand at (b) 50 °C. Each tabulated value in the contact maps is calculated by averaging over the actual number of contacts of that type obtained from the four model systems and dividing by the number of that type that would exist if the same total number of contacts were assigned randomly on the basis of the fraction of each monomer type present. The corresponding standard errors of the two contact maps are shown in (c) and (d).

Figure 13. Snapshots of METHOCEL A 10-mers in water at 4.5 wt % (left) and in acetone at 6.4 wt % (middle) at the end of a 25 ns simulation and oligomer−oligomer RDFs of METHOCEL A 10-mers in water and acetone (right).

one of the monomers is less than 0.4 nm from any atom on the other. For the “10-mer 4 wt %”, “20-mer 3 wt %”, “20-mer 6 wt %”, and “40-mer 4 wt %” systems studied here, three repeat simulations with random initial configurations are set up and simulated for 35 ns. A 200 ps production run is then conducted for each simulation, and five frames separated by 50 ps intervals are used to compute the average number of contacts for each of the 64 types of monomer−monomer pairs, where the monomer “types” are distinguished by the eight methyl substitution patterns ranging from no methyl substituents to trisubstituted MC. Given the fractions of each type of monomer present and the total number of contacts, we compute the numbers of each pair of monomer types that would be expected if the contacts between them were random, and we use this to normalize the raw counts of actual contact numbers. A value greater than unity then implies that the

contact occurs more often than would be expected if monomers form contacts randomly with other monomers. We then average these normalized counts over five frames, divide by the number of monomers in the system, and average the results over three simulations for each oligomer system to produce the contact map for each of the four systems studied at two different temperatures, giving the eight maps shown in Figure SI5a,b in the Supporting Information. At each temperature, we observe that the majority of the contact values (∼90%) are within the standard error among the four systems. We therefore further average the values in the four contact maps at the same temperature to produce contact maps averaged over the four oligomer model systems at 25 and 50 °C, which are shown in Figure 12a,b respectively. The corresponding standard errors, presented in Figure 12c,d, are generally below 0.1, indicating that the generated contact maps are largely insensitive to the N



■

Article

CONCLUSIONS We used two GROMOS force fields to simulate single cellulose and methylcellulose oligomers with lengths ranging from dimers to 40-mers in water. Despite some differences in conformational distributions at the dimer level, the two force fields give similar predictions for the stiffness of longer chains and the numbers of hydrogen bonds for various methyl substitution levels and locations. We then used the newer force field of the two, GROMOS 56Acarbo, to simulate multiple methylcellulose homo-oligomers in water at low concentrations (around 4.5 wt %) at 25 and 90 °C. We found that the three substitution positions have different hydrogen-bond-forming preferences. The O6 position forms the most hydrogen bonds, followed by O2 and O3. Dimethylated and trimethylated MCs form loosely aggregated structures at low temperature (25 °C), and most MCs form precipitate-like structures at high temperature (90 °C). O6−O6 is the most abundant intermolecular oligomer−oligomer hydrogen bond formed in many MC systems and, with the presence of the other intermolecular oligomer−oligomer hydrogen bonds, is responsible for the different solvation behaviors of methylcellulose molecules with similar degrees of substitution. We created model random METHOCEL A copolymers with substitution patterns mimicking those of the commercial system and studied their gelation at both room temperature and 50 °C, which is near the experimentally determined gelation temperature. We found that short oligomer chains suffice to represent the temperature dependence of RDF peak heights and the hydrogen bonding of chains up to 4 times longer. The results for our model oligomer systems also suggest that hydrophobic interactions induce methylcellulose gelation, which supports one of the proposed mechanisms reported in the literature. However, even with state-of-the-art computational power, a 22 nm cubic simulation box, and a 40-mer model oligomer chain, we cannot simulate the fibril formation seen in experimental methylcellulose gels. Finally, we applied the same force field (GROMOS 56Acarbo) to simulate these oligomers in acetone and observed MC oligomer dispersal rather than the aggregation seen in water, as expected.

randomly generated sequence based on the statistics of the commercial METHOCEL A products. A few contact-map values involving cellulose, 3-MC, and 3,6-MC, however, have relatively large standard errors, likely due to the low probabilities of their presence in the METHOCEL A product (mole fractions ≤ 0.05). We use the “hot spots” (i.e., high normalized contact numbers, marked in red or purple) on the contact maps to suggest which monomer pairs are most responsible for inducing gelation at elevated temperature. The normalized count of the “2,3,6-MC−2,3,6-MC” pair is the highest among all pairs at both temperatures, and the count at elevated temperature is larger than that at room temperature. This suggests that hydrophobic interactions between the trisubstituted monomers are responsible for the formation of aggregates at temperatures above the experimentally determined gelation temperature. The intensity of the “Un−Un” (cellulose−cellulose) pair also increases at elevated temperature. This observation correlates well with the increased number of intermolecular hydrogen bonds at elevated temperature shown in Table 7. However, the “Un−Un” pair intensities do not exceed the average pair contact intensity (e.g., intensity < 1.0), suggesting that hydrogen bonding is not the main driving force for the formation of aggregates in these methylcellulose oligomer systems. Our observations from these contact maps agree with the methylcellulose gelation mechanism proposed by Kato et al.11 and Li et al.,12 in which trisubstituted methylcellulose units act as hydrophobic junctions. In addition to simulating METHOCEL A model oligomer chains in water, we also perform simulations in acetone. Figure 13, which shows final snapshots of 10-mer model oligomers at concentrations of 4.5 and 6.4 wt %, respectively, in water and acetone, clearly differentiates the solvation behavior in these two solvents. The oligomer−oligomer RDF indeed shows a drastic decrease in peak height in acetone relative to water, indicating that no oligomer cluster formation occurs in the former. To our knowledge, this is the first attempt to simulate multiple methylcellulose oligomer chains in an organic solvent. This is of practical interest because METHOCEL soliddispersion drug preparations are prepared in organic solvents such as acetone and acetic acid. Therefore, the creation of a set of model METHOCEL A chains that can be simulated in different solvents will allow us to explore the interactions between oligomers and drugs in those solvents at an atomistic level. Our study is the first to model methylcellulose gelation behavior using atomistic-level MD simulations. We have shown that all four model METHOCEL A systems form closer-packed clusters at an elevated temperature around the experimentally determined gelation temperature than they do at lower (room) temperature. While the currently attainable chain lengths, box sizes, and time scales are not sufficient for atomistic simulations of the fibril gelation structure observed experimentally in methylcellulose gels, even with short-length model oligomers (10-mers) we are able to represent the solvation behavior of the METHOCEL A chemistry through RDFs and hydrogen-bond counts and to probe the gel-formation mechanism through contact maps. Therefore, despite the limitations of our method, it is possible to gain insight into the interaction between these polymers and small molecules such as drugs or peptides.

■

ASSOCIATED CONTENT

S Supporting Information *

Additional data for single-chain Rg values, RDF curves, percentages of intermolecular hydrogen bonds, and contact maps for each of the four model METHOCEL A systems with the corresponding standard errors. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address §

(I.S.D.) Department of Chemical Engineering, Indian Institute of Technology Kanpur, Kanpur, Uttar Pradesh 208016, India. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported in part through computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor, and in part through the National Science Foundation (Grant OCIO



Article

(19) Wang, D.; Á mundadóttir, M. L.; van Gunsteren, W. F.; Hünenberger, P. H. Intramolecular Hydrogen-Bonding in Aqueous Carbohydrates as a Cause or Consequence of Conformational Preferences: A Molecular Dynamics Study of Cellobiose Stereoisomers. Eur. Biophys. J. 2013, 42, 521−537. (20) Yu, H.; Amann, M.; Hansson, T.; Köhler, J.; Wich, G.; van Gunsteren, W. F. Effect of Methylation on the Stability and Solvation Free Energy of Amylose and Cellulose Fragments: A Molecular Dynamics Study. Carbohydr. Res. 2004, 339, 1697−1709. (21) Winger, M.; Christen, M.; van Gunsteren, W. F. On the Conformational Properties of Amylose and Cellulose Oligomers in Solution. Int. J. Carbohydr. Chem. 2009, 307695. (22) Shen, T.; Langan, P.; French, A. D.; Johnson, G. P.; Gnanakaran, S. Conformational Flexibility of Soluble Cellulose Oligomers: Chain Length and Temperature Dependence. J. Am. Chem. Soc. 2009, 131, 14786−14794. (23) Casay, G. A.; George, A.; Hadjichristidis, N.; Lindner, J. S.; Mays, J. W.; Peiffer, D. G.; Wilson, W. W. Dilute Solution Properties, Chain Stiffness, and Liquid Crystalline Properties of Cellulose Propionate. J. Polym. Sci., Part B: Polym. Phys. 1995, 33, 1537−1544. (24) Queyroy, S.; Müller-Plathe, F.; Brown, D. Molecular Dynamics Simulations of Cellulose Oligomers: Conformational Analysis. Macromol. Theory Simul. 2004, 13, 427−440. (25) Srinivas, G.; Cheng, X.; Smith, J. C. A Solvent-Free Coarse Grain Model for Crystalline and Amorphous Cellulose Fibrils. J. Chem. Theory Comput. 2011, 7, 2539−2548. (26) Hansen, H. S.; Hünenberger, P. H. A Reoptimized GROMOS Force Field for Hexopyranose-Based Carbohydrates Accounting for the Relative Free Energies of Ring Conformers, Anomers, Epimers, Hydroxymethyl Rotamers, and Glycosidic Linkage Conformers. J. Comput. Chem. 2011, 32, 998−1032. (27) Lins, R. D.; Hünenberger, P. H. A New GROMOS Force Field for Hexopyranose-Based Carbohydrates. J. Comput. Chem. 2005, 26, 1400−1412. (28) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43−56. (29) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (30) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718. (31) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. Intermol. Forces 1981, 14, 331−342. (32) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126, No. 014101. (33) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (34) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695−1697. (35) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (36) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (37) Nosé, S.; Klein, M. L. Constant Pressure Molecular Dynamics for Molecular Systems. Mol. Phys. 1983, 50, 1055−1076. (38) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. (39) Kräutler, V.; Müller, M.; Hünenberger, P. H. Conformation, Dynamics, Solvation and Relative Stabilities of Selected B-Hexopyranoses in Water: A Molecular Dynamics Study with the GROMOS 45A4 Force Field. Carbohydr. Res. 2007, 342, 2097−2124. (40) Nishida, Y.; Ohrui, H.; Meguro, H. 1H-NMR Studies of (6r)and (6s)-Deuterated d-Hexoses: Assignment of the Preferred

1053575) Extreme Science and Extreme Science and Engineering Discovery Environment (XSEDE) resources provided by the XSEDE Science Gateways Program (XSEDE Grant TGCHE140009). The work was funded by The Dow Chemical Company. We are also grateful for very useful discussions with William “Trey” Porter of Dow Chemical Co. and Prateek K. Jha.

■

REFERENCES

(1) U.S. Food and Drug Administration. GRAS Substances (SCOGS) Database. http://www.fda.gov/Food/ IngredientsPackagingLabeling/GRAS/SCOGS/default.htm (accessed Sept 4, 2013). (2) METHOCEL Cellulose Ethers Technical Handbook; The Dow Chemical Company: Midland, MI, 2002. (3) Siepmann, J.; Peppas, N. A. Modeling of Drug Release from Delivery Systems Based on Hydroxypropyl Methylcellulose (HPMC). Adv. Drug Delivery Rev. 2001, 48, 139−157. (4) Ilevbare, G. A.; Liu, H.; Edgar, K. J.; Taylor, L. S. Inhibition of Solution Crystal Growth of Ritonavir by Cellulose PolymersFactors Influencing Polymer Effectiveness. CrystEngComm 2012, 14, 6503− 6514. (5) Keary, C. M. Characterization of METHOCEL Cellulose Ethers by Aqueous SEC with Multiple Detectors. Carbohydr. Polym. 2001, 45, 293−303. (6) Klemm, D.; Heublein, B.; Fink, H.-P.; Bohn, A. Cellulose: Fascinating Biopolymer and Sustainable Raw Material. Angew. Chem., Int. Ed. 2005, 44, 3358−3393. (7) Beer, M. U.; Wood, P. J.; Weisz, J. A Simple and Rapid Method for Evaluation of Mark−Houwink−Sakurada Constants of Linear Random Coil Polysaccharides Using Molecular Weight and Intrinsic Viscosity Determined by High Performance Size Exclusion Chromatography: Application to Guar Galactomannan. Carbohydr. Polym. 1999, 39, 377−380. (8) Sarkar, N. Thermal Gelation Properties of Methyl and Hydroxypropyl Methylcellulose. J. Appl. Polym. Sci. 1979, 24, 1073− 1087. (9) Nilsson, S.; Sundelöf, L.-O.; Porsch, B. On the Characterization Principles of Some Technically Important Water Soluble Non-Ionic Cellulose Derivatives. Carbohydr. Polym. 1995, 28, 265−275. (10) Kobayashi, K.; Huang, C.; Lodge, T. P. Thermoreversible Gelation of Aqueous Methylcellulose Solutions. Macromolecules 1999, 32, 7070−7077. (11) Kato, T.; Yokoyama, M.; Takahashi, A. Melting Temperatures of Thermally Reversible Gels IV. Methyl Cellulose−Water Gels. Colloid Polym. Sci. 1978, 256, 15−21. (12) Li, L.; Thangamathesvaran, P. M.; Yue, C. Y.; Tam, K. C.; Hu, X.; Lam, Y. C. Gel Network Structure of Methylcellulose in Water. Langmuir 2001, 17, 8062−8068. (13) Haque, A.; Morris, E. R. Thermogelation of Methylcellulose. Part I: Molecular Structures and Processes. Carbohydr. Polym. 1993, 22, 161−173. (14) Fairclough, J. P. A.; Yu, H.; Kelly, O.; Ryan, A. J.; Sammler, R. L.; Radler, M. Interplay between Gelation and Phase Separation in Aqueous Solutions of Methylcellulose and Hydroxypropylmethylcellulose. Langmuir 2012, 28, 10551−10557. (15) Bodvik, R.; Dedinaite, A.; Karlson, L.; Bergströ m, M.; Bäverbäck, P.; Pedersen, J. S.; Edwards, K.; Karlsson, G.; Varga, I.; Claesson, P. M. Aggregation and Network Formation of Aqueous Methylcellulose and Hydroxypropylmethylcellulose Solutions. Colloids Surf., A 2010, 354, 162−171. (16) Lott, J. R.; McAllister, J. W.; Arvidson, S. A.; Bates, F. S.; Lodge, T. P. Fibrillar Structure of Methylcellulose Hydrogels. Biomacromolecules 2013, 14, 2484−2488. (17) Matthews, J. F.; Himmel, M. E.; Brady, J. W. Simulations of the Structure of Cellulose. ACS Symp. Ser. 2010, 1052, 17−53. (18) Bergenstråhle-Wohlert, M.; Brady, J. W. Overview of Computer Modeling of Cellulose. Methods Mol. Biol. 2012, 908, 11−22. P



Article

Rotamers about C5−C6 Bond of D-Glucose and D-Galactose Derivatives in Solutions. Tetrahedron Lett. 1984, 25, 1575−1578. (41) Ajisaka, K.; Fujimoto, H.; Nishida, H. Enzymic Synthesis of Disaccharides by Use of the Reversed Hydrolysis Activity of β-Dgalactosidases. Carbohydr. Res. 1988, 180, 35−42. (42) Stenutz, R.; Carmichael, I.; Widmalm, G.; Serianni, A. S. Hydroxymethyl Group Conformation in Saccharides: Structural Dependencies of 2JHH, 3JHH, and 1JCH Spin−Spin Coupling Constants. J. Org. Chem. 2002, 67, 949−958. (43) Patel, T. R.; Morris, G. A.; de la Torre, J. G.; Ortega, A.; Mischnick, P.; Harding, S. E. Molecular Flexibility of Methylcelluloses of Differing Degree of Substitution by Combined Sedimentation and Viscosity Analysis. Macromol. Biosci. 2008, 8, 1108−1115. (44) Hiemenz, P. C.; Lodge, T. P. Polymer Chemistry, 2nd ed.; CRC Press: Boca Raton, FL, 2007. (45) van der Spoel, D.; van Maaren, P. J.; Larsson, P.; Tîmneanu, N. Thermodynamics of Hydrogen Bonding in Hydrophilic and Hydrophobic Media. J. Phys. Chem. B 2006, 110, 4393−4398. (46) van Gunsteren, W. F.; Billeter, S. R.; Eising, A. A.; Hünenberger, P. H.; Krüger, P.; Mark, A. E.; Scott, W. R. P; Tironi, I. G. Biomolecular Simulation: The GROMOS96 Manual and User Guide; Verlag der Fachvereine Hochschulverlag AG an der ETH Zürich: Zürich, Switzerland, 1996. (47) Srinivas, G.; Klein, M. L. Molecular Dynamics Simulations of Self-Assembly and Nanotube Formation by Amphiphilic Molecules in Aqueous Solution: A Coarse-Grain Approach. Nanotechnology 2007, 18, No. 205703.

Q


Analysis of Solvation and Gelation Behavior of Methylcellulose Using

Recommend Documents