Article pubs.acs.org/Macromolecules
A Systematic Coarse-Grained Model for Methylcellulose Polymers: Spontaneous Ring Formation at Elevated Temperature Wenjun Huang, Rahul Ramesh, Prateek K. Jha, and Ronald G. Larson* Department of Chemical Engineering, University of Michigan, Ann Arbor, Michigan 48109-2136, United States S Supporting Information *
ABSTRACT: We develop a systematic coarse-grained (CG) model for methylcellulose polymers, including random copolymers with compositions representative of modeling commercial METHOCEL A polymer, using one CG bead per monomer. We parametrize our CG model using the RDFs from atomistic simulations of short methylcellulose oligomers, extrapolating the results to long chains. Using a LJ 9−6 potential, the CG model captures the effect of monomer substitution type and temperature observed in detailed atomistic simulations. We use dissociation free energy to validate our CG model against the atomistic model. We then use this CG model to simulate single chains up to 1000 monomers long, and we calculate persistence lengths for a selection of homogeneous and heterogeneous methylcellulose chains, which show good agreement with experimental results. Interestingly, simulations of 600-mer heterogeneous chains show a collapse transition at 50 °C and form a stable ring structure with outer diameter around 14 nm. This structure appears to be a precursor to fibril structure reported in a recent study of methylcellulose gels [Biomacromolecules 2013, 14, 2484].
■
INTRODUCTION Commercial methylcellulose (MC) products developed by the Dow Chemical Company, marketed under the brand name METHOCEL A, are widely used in agricultural, ceramic processing, construction, and pharmaceutical industries. MC is categorized as safe to use as a food additive by the U.S. Food and Drug Administration.2 Besides its commercial value, MC is of considerable scientific interest as a self-attractive semiflexible water-soluble random copolymer that forms a gel at elevated temperatures and has been studied both experimentally and theoretically.3,4 To form the MC monomer, up to three reactive hydroxyl groups (−OH) on the natural cellulose monomer are substituted with hydrophobic methyl groups (−CH3). Methylcellulose monomer, as a result of this substitution, can have a range of hydrophobicities and degrees of substitution (DS), defined as the moles of substituents (i.e., hydrophobic methyl groups) per mole of MC monomer.5 For convenience, we will here take the term “methylcellulose (MC) monomer” to include the cellulose monomer as a limiting case with DS = 0. We will also use the term “cellulosic” to encompass both unmethylated and methylated monomers and polymers throughout the paper. Commercial METHOCEL A product has two important properties: it is soluble in water and can form a thermoreversible gel at elevated temperatures. Fully substituted MC (i.e., DS = 3) is insoluble in water because of its strong hydrophobic interactions. Cellulose (DS = 0) is also insoluble in water as the presence of multiple hydroxyl groups on a monomer gives rise to intramolecular hydrogen bonds and results in cellulose crystallization. However, MC with DS around 2 is water-soluble as a result of the balanced effect of breaking the tendency of © 2016 American Chemical Society
cellulose to crystallize via incorporation of hydrophobic methyl groups while retaining enough hydrophilic hydroxyl groups to sustain enough hydrogen bonding with water.6 The gel morphology that forms in MC gels at elevated temperatures has recently been identified by Lott et al.1 using cryo-TEM as a network of fibril structures with a uniform diameter of around 14 ± 2 nm above 55 °C. The fibril formation mechanism, however, is still unclear. Previous atomistic MD simulations in our group7 of MC oligomers suggested that the hydrophobic interaction is the major driving force for oligomers to form clusters, in line with the earlier theory proposed by Kato et al.8 Nevertheless, these earlier simulations and theoretical works could not provide a satisfactory answer regarding why the diameter of the fibril stops increasing beyond 14 nm. Collapse transitions of self-attractive semiflexible polymers have been well-documented in experimental studies.9,10 Examples of self-attractive semiflexible polymers include biopolymers such as F-actin and DNA and synthetic polymers such as Kevlar and Zylon used in high-performance fibers. Many simulation studies11−14 have been carried out over the past decade to characterize the collapse transitions of these selfattractive semiflexible polymers under different solvent conditions. A recent systematic simulation study by Kong and Larson,13 for example, detailed various collapsed states and collapse paths exhibited by semiflexible polymers with different persistence lengths and attractive interaction strengths. Other important simulation studies have shown similar collapsed Received: October 31, 2015 Revised: January 6, 2016 Published: February 11, 2016 1490
DOI: 10.1021/acs.macromol.5b02373 Macromolecules 2016, 49, 1490−1503
Article
Macromolecules states.15−18 These transitions often involve long-lived intermediate states such as “hairpin structures” and final compact states including toruses, condensed globules, and folded bundles. These previous simulation studies of semiflexible polymers utilized a top-down approach, involving generic bead−spring models that were not targeted to any specific polymer chemistry. Our paper describes a bottom-up approach, with the force field designed to represent the specific physicochemical properties of a particular semiflexible polymer. While the CG force field is specialized to the polymer of interest, our method of developing it is general, and the structures we report, especially rings or toruses, have shown up in simulations of generic models, indicating their likely occurrence in collapse transitions of other semiflexible polymers. Several computational simulation techniques have been used to study the solvation and gelation behavior of cellulose and its derivatives, including methylcellulose. Unmethylated cellulose oligomers and cellulose melts have also been simulated.19−22 The solvation behavior of multiple cellulosic oligomers in water and organic solvents has been explored,7 and interactions between cellulosic oligomers and small drug molecules have also been investigated. 23 Even though these atomistic simulations provide important information regarding local chain orientation and intermolecular interaction at the nanometer scale for cellulosic polymers, the time and length scales that could be simulated were too small to reveal the structure and dynamics of gel formation. Thus, there is a need for a coarse-grained (CG) force field to understand the structure and dynamics of these cellulosic polymers at larger length and time scales, while retaining sufficient chemicalspecific information to distinguish the interactions of monomers with differing degrees of methyl substitution. A number of CG force fields have been developed for crystalline or amorphous cellulose structures. The MARTINI and the M3B force fields24,25 both adopted three-site CG models (i.e., wherein each cellulose monomer is represented by three CG beads). Srinivas et al.26 adopted a one-site CG model and used it to demonstrate the transition of cellulose from a crystalline fibril to an amorphous state. In this work, we present an implicit-solvent CG force field parametrized for MC, including all eight monomer substitution types ranging from DS = 0 to 3, at both room temperature (25 °C) and an elevated temperature (50 °C). The force field is developed by matching the monomer−monomer radial distribution functions (RDFs) obtained using the CG model to corresponding results obtained using atomistic force field for short homo-oligomer chains and extrapolating the resulting force field parameters to longer chains. Using Brownian dynamics (BD) simulations, we show that the CG force field can match the oligomers’ dissociation free energies obtained from corresponding atomistic simulations qualitatively. The force field can also distinguish the dependence of polymer selfassociation on the temperature and on monomer substitution type effect among different MC monomers. We apply the force field to both homogeneous and heterogeneous MC chains and obtain chain persistence lengths that are in agreement with experimental results. Our simulations reveal the formation of a ring structure for highly substituted (DS ≥ 2) long chains, which we believe is the precursor to the methylcellulose fibril structure.1 We also calibrate the time scale of the CG BD simulations by comparison of chain relaxation rates to those from polymer theory. This work demonstrates that parameters
obtained from systematic coarse graining can be used to study macroscopic structural properties and aggregation behavior of realistic polymer chains that are much longer than their persistence length. Our approach can be adapted to polymers used in drug delivery applications where the relevant size of molecules is beyond the reach of atomistic MD techniques. In addition, due to the use of an analytical LJ potential, the number of parameters in our model is limited. For homopolymer chains the parameters can be reduced to a dimensionless interaction strength and a ratio of persistence length to effective polymer diameter. Thus, the results of our model can be fit to a general phase diagram describing the collapsed states of semiflexible self-attracting polymers. The paper is organized as follows. We first discuss our coarse-graining strategy and the choice of parameters. We then present persistence length estimations of both homogeneous and heterogeneous MC chains based on simulated radius of gyration (Rg) values obtained using the CG parameters. Finally, we show results of long chain simulations of MC polymers that demonstrates the formation of ring structures.
■
METHODS Atomistic Simulation. The atomistic molecular dynamics (MD) simulation technique used in this work has been described in detail in an earlier study,7 and here we will mention only a few details briefly. The simulations were carried out in GROMACS,27−29 version 4.6.5, using the GROMOS 56Acarbo force field.30 Materials Studio, version 7.0.1 (Accelrys Inc.), was used to build MC monomers with all possible combinations of methyl substituents. Production runs were carried out using an NPT ensemble with a Nosé−Hoover thermostat31,32 and a Parrinello−Rahman barostat33,34 with a time step of 1 fs. Equilibrium atomistic MD simulations were used for CG force field parametrization. For 5-, 10-, 15-, and 20-mer systems, cubic simulation boxes of 10 nm length were used, while for 25- and 30-mer systems a 15 nm cubic simulation box was used. The polymer concentration was kept at 10 wt % unless otherwise noted. Periodic boundary conditions were employed on all three directions. The simulation box was solvated with SPC water molecules which are compatible with the GROMOS force field. Production runs were carried out for at least 20 ns, and the last 5 ns was used to obtain radial distribution functions (RDFs). Simulations were conducted at two different temperatures (25 and 50 °C). “Pulling” simulations were used to extract dissociation free energy between two chains, where two 10-mer chains were placed in parallel in a 5 nm × 8 nm × 13 nm simulation box. The z component of the center-of-mass (COM) distance (DCOM_Z) between the two chains was calculated and used as the pulling coordinate. A spring constant of 1000 kJ−1 nm−2 and a pulling rate of 0.001 nm ps−1 were used to generate a pulling trajectory between DCOM_Z = 0.5 nm and DCOM_Z = 4.2 nm. In this range of DCOM_Z values, a series of configurations were chosen at every 0.1 nm increment, resulting in a total of 39 configurations. Each configuration was used to start a 40 ns MD simulation for umbrella sampling, where a harmonic force constant of 1300 kJ mol−1 nm−2 was applied throughout the simulation. A weighted histogram analysis method (WHAM)35 implemented using the GROMACS utility g_wham36 was used to generate three potential of mean force (PMF) diagram. Coarse-Grained (CG) Simulations. We chose a Brownian dynamics (BD) simulation technique to perform the CG simulations. The simple BD method used here, although it 1491
DOI: 10.1021/acs.macromol.5b02373 Macromolecules 2016, 49, 1490−1503
Article
Macromolecules
Figure 1. Schematics of the methylcellulose coarse-grained model. Each methylcellulose monomer (DS ranging from 0 to 3) is represented by one bead centered at the monomer center-of-mass (COM). The beads are connected via hard harmonic springs.
Figure 2. (a) Intramolecular atomistic monomer RDFs obtained from atomistic simulations of 10-mer single chain homogeneous methylcellulose. The RDF is an average of RDFs obtained from eight single homopolymer chain systems. Each homopolymer consists of a chain of 10 monomers of one of the eight methylcellulose monomer substitution types. (b), (c), and (d) are zoomed-in view of the second, third, and fourth peak of those shown in (a) at 1.0, 1.5, and 2.0 nm, respectively.
simulation package37 (ver. Feb 2014) in an NVE ensemble. Simulations were set up with dimensionless LJ (Lennard-Jones) units, with the Boltzmann constant (kb) and the three fundamental units (scales) defined as the mass (m), distance (σ), and energy (ε). The dimensionless values of particle mass, bond length, and temperature are taken equal to unity. The conversion factors between these fundamental dimensionless
neglects long-range hydrodynamic interactions, is a suitable tool for this study because we are primarily interested in the interaction between hydrophobic polymers rather than interaction between polymer and solvent, and the effect of solvent on polymer−polymer interactions can be captured implicitly through the CG effective polymer−polymer interactions. BD simulations were performed using LAMMPS 1492
DOI: 10.1021/acs.macromol.5b02373 Macromolecules 2016, 49, 1490−1503
Article
Macromolecules
height (g(r) value) of the first, second, and third peak of the same atomistic RDF. We found that the intramolecular atomistic RDFs are similar among all 10-mer homo-MC chains with different monomer substitution types, which is expected because the contour length of 10-mer chains (∼5 nm) is well below the persistence length of MC (∼11 nm),40 and therefore these stiff chains do not show the effect of substitution. On the basis of the information from the intramolecular atomistic monomer RDF, we decided to average the RDFs obtained for all eight homo-MC chains and to use a single set of bonded parameters for all MCs. Note that when methylcellulose monomer is substituted at the 3-position (i.e., 3-MC), the intrachain hydrogen bonding network is disrupted and therefore the chain is more flexible.19 Yet this effect is not predominant in a MC chain that is shorter than 40 monomers long. Therefore, we choose to capture this effect by tuning the nonbonded interaction, rather than bonded interaction. A summary of all bonded parameters is tabulated in Table 1.
quantities and their dimensional counterparts have been determined and will be discussed later in the text. Using a Langevin thermostat,38 the temperature was maintained at 1 dimensionless temperature unit (298 K). Note here we kept the dimensionless temperature parameter constant regardless of the “effective temperature” at which we are running the simulation. The “effective temperature” is implicitly captured through the intermolecular interaction parameters parametrized from atomistic simulations conducted at different temperatures. The damping coefficient was set to 10 and a typical time step of 0.001τ, where τ = (mσ2/ϵ)1/2 being the time scale. A typical simulation was carried out for at least 107 steps depending on the length of the chain, and the final one-third of the data were used to obtain the averaged Rg values for the polymer chain unless the Rg values converged to a specific value earlier than this, which occurred in the case of chain collapse. Because solvent molecules were implicitly represented, the size of the box was arbitrarily chosen to be 1.5 times the contour length of the longest polymer chain in the box. CG “pulling” simulations were set up in a similar fashion to the atomistic counterparts. The pulling distance was from 1.0σ to 7.0σ; the harmonic spring force constant for each histogram was set to 7.0ε, and a total of 30 simulation windows with 0.2σ spacing were used to generate PMF. A WHAM calculation code (V2.0.9) developed by Grossfield lab was used in the PMF calculation.39 Coarse-Grained (CG) Model. The CG polymer chains were modeled using beads and stiff springs, with each bead located at the center-of-mass (COM) of a MC monomer, as shown in Figure 1. We included both bonded and nonbonded interactions in our bead/stiff-spring model. Because each bead represents one charge-neutral MC monomer, we did not include any explicit electrostatic interaction terms. The complete CG interaction potential was expressed by the equation UCG = Ubond + Uangle + Udihedral + Unonbonded
Table 1. Summary of the Intramolecular Parameters in CG Modela Bonded Parameters dimensionless units l0 kbond θ0 kangle n, d kdihedral
1 1000 165 30 1 2
dimensional units
σ ε/σ2 deg ε/rad2
0.515 2478.28 2.88 74.35
ε 4.96 Nonbonded Parameters dimensionless units
εii,long‑range σii,long‑range rc_ii,long‑range
(1)
nm kJ/(mol nm2) rad kJ/(mol rad2)
0.2 0.874 4
ε σ σ
kJ/mol dimensional units 0.5 0.450 2.06
kJ/mol nm nm
εii, σii, and rc_ii are the input values to the Lennard-Jones 9−6 potential (eq 5). a
The bonded interactions included harmonic bond, angle, and dihedral interactions, which were applied to any two, three, and four consecutive beads on a chain, respectively (eqs 2−4). Ubond =
1 Kb(l − l0)2 2
(2)
Uangle =
1 K θ(θ − θ0)2 2
(3)
Udihedral = Kφ[1 + d cos(nφ)]
We selected the conversion factors between the three fundamental dimensionless units (mass m, distance σ, and energy ε) and their dimensional units counterparts. The unit mass was chosen to be 188 Da, which is the averaged molecular weight of all eight MC monomers. The unit length was set to be 0.515 nm, which is the simulated average COM separation of monomers from the atomistic RDF. The unit energy was chosen to be 2.478 kJ mol−1, corresponding to 1kBT at 298 K. The converted bonded parameters in real units are also tabulated in Table 1. All other dimensionless units can be expressed as combinations of these fundamental dimensionless units, as detailed in the LAMMPS manual, and therefore can be also converted to the dimensional units. The nonbonded interactions took the form of a truncated and shifted Lennard-Jones (LJ) 9−6 potential (eq 5). The choice of the LJ 9−6 potential instead of the more commonly adopted LJ 12−6 potential reflects the flat geometry of the methylcellulose monomer. The potential was shifted upward so that the value of the potential function goes to zero at a cutoff distance (rc). We introduced an additional weak intramolecular nonbonded interaction between any four and five consecutive beads on one chain to fine-tune the CG RDF so that the remaining peak positions (fourth and fifth peaks) on the intramolecular atomistic monomer RDFs were matched with
(4)
Here l0 and θ0 are the equilibrium bond length and angle, and Kb and Kθ are the corresponding bond and angle force constants, respectively. In the dihedral expression, d and n are the phase constants, and Kφ is the dihedral force constant. These parameters were determined from a single 10-mer chain atomistic simulation by mapping the intramolecular CG bead− bead radial distribution function (RDF) onto the corresponding atomistic single chain intramolecular monomer COMmonomer COM RDF, referred to hereafter as the “intramolecular atomistic monomer RDF”. A typical fit of the CG RDF to the intramolecular atomistic monomer RDF is shown in Figure 2. The equilibrium bond length (l0) and angle (θ0) were determined by matching the peak position (r value) of the first (∼0.5 nm) and the second peak (∼1.0 nm) of the intramolecular atomistic RDF. The bond, angle, and dihedral constants (Kb, Kθ, and Kφ) were determined by matching the 1493
DOI: 10.1021/acs.macromol.5b02373 Macromolecules 2016, 49, 1490−1503
Article
Macromolecules
Figure 3. (a) Intermolecular atomistic monomer RDFs for 10 wt % 2,6-MC with different chain lengths at 25 °C. To reduce statistical noise, we averaged RDFs of 5- and 10-mers (red), 15- and 20-mers (green), and 25- and 30-mers (blue). (b) Same as (a) except these RDFs were generated at 50 °C.
and extrapolate these to longer chains. This can only be achieved through the use of an analytical potential function. (3) Commercial MC products contain all eight MC monomer substitution types. Using an analytical potential function for each monomer substitution type is particularly convenient for simulating heterogeneous MC chains, where a geometric mixing rule is used to calculate the εij and σij between different pairs (eq 6). If, however, tabulated potentials were to be employed, the 24 cross terms between all eight MC monomer substitution combination types would need to be computed at both short and long chain lengths, which would require a tremendous modeling effort.
CG RDFs. Because these RDFs were similar among all MC monomer substitution types, these weak nonbonded 1−4 and 1−5 interactions were kept the same for all MCs and are tabulated in Table 1. ⎧ ⎡ 9 ⎛ σ ⎞6 ⎤ ⎛ σ ⎞ 9 ⎛ σ ⎞6 ⎛ σ ⎞ ⎪ ⎪ 4εii⎢⎜ ii ⎟ − ⎜ ii ⎟ − ⎜ ii ⎟ + ⎜ ii ⎟ ⎥ r < rc ⎝r ⎠ Unonbonded(r ) = ⎨ ⎢⎣⎝ r ⎠ ⎝ rc ⎠ ⎝ rc ⎠ ⎥⎦ ⎪ ⎪0 r ≥ rc ⎩
(5)
The intermolecular nonbonded interaction parameters were fitted from intermolecular monomer COM−monomer COM RDFs, which are referred as “intermolecular atomistic monomer RDFs” in the following. These were generated from atomistic simulations with 10 wt % polymer loading at various chain lengths. These intermolecular RDFs were generated in a similar way as the intramolecular RDFs, except that the contributions from the five neighboring monomers of any given monomer were omitted. A typical intermolecular atomistic monomer RDF had several closely spaced peaks at short distance (see Figure 3). Achieving a good fit of the entire RDF of this kind would require a tabulated potential, as demonstrated by Srinivas et al.26 We, however, decided to use an analytical form, namely the LJ 9−6 potential, even though this choice allowed us to achieve good fits of only the first main peak, which occurs roughly at r = 0.6 nm in all intermolecular atomistic monomer RDFs. Our choice was justified by the following reasons: (1) The first peak of the intermolecular atomistic monomer RDF reflects the equilibrium COM distance between two monomers on different chains, which correlates with the intermolecular interaction strength between these two monomers, and is captured by the 9−6 LJ potential in our CG force field. (2) Our goal is to simulate MC chains with realistic chain lengths (>400 monomers), which are well beyond the longest chain length we can afford to simulate at an atomistic scale (30 monomers). As a result, obtaining tabulated potentials for such long chains would be extremely challenging due to the lack of reference atomistic simulations. Our strategy is to use an analytical expression to obtain fit parameters for short chains
σij =
σiiσjj
εij =
εiiεjj
(6)
(4) Our use of analytical potentials in a CG model can be adapted to model other related and important experimental polymers, including hydroxypropyl methylcellulose (HPMC) and its derivatives. These polymers are often used as drug carriers, and the sizes of the polymer−drug complexes are usually around a few hundred nanometers. We used a modified iterative Boltzmann inversion (IBI) scheme to obtain the intermolecular nonbonded interaction parameters. The standard IBI scheme (eq 7) has been used in many studies to achieve very good fitting between CG and atomistic RDFs:26,41−44 ⎛ g (r ) ⎞ Vi + 1(r ) = Vi (r ) + αkBT ln⎜⎜ i ⎟⎟ ⎝ g (r ) ⎠
(7)
In this equation, Vi(r) and Vi+1(r) are the tabulated potentials in the ith and (i + 1)th iteration, respectively; gi(r) and g(r) are the RDFs corresponding to the potential at the ith iteration and target RDF, respectively. α is the damping factor which is arbitrarily chosen and decreases after each iteration. We modified the standard IBI scheme and used an analytical potential (eqs 8 and 9) to substitute for the tabulated potential form of the ith iteration. Ui(r ) = εif (r )
(8)
where 1494
DOI: 10.1021/acs.macromol.5b02373 Macromolecules 2016, 49, 1490−1503
Article
Macromolecules
Figure 4. Fits of the intermolecular CG RDFs to intermolecular atomistic monomer RDFs for four homomethylcellulose systems. The atomistic systems are cellulose, 3-MC, 2,6-MC, and 2,3,6-MC. The atomistic RDFs were obtained by computing the intermolecular monomer COM− monomer COM RDF for each 10 wt % 15-mer homo-oligomer system. The CG RDFs were computed by the same method were fit to match the position (rp) and height (g(rp)) of the first peak in atomistic RDFs.
⎛ (r − r )2 ⎞ p ⎟ w(r ) = exp⎜⎜ − ⎟ 2 2 s ⎝ ⎠ d
⎧ ⎡ 9 ⎛ σ ⎞6 ⎤ ⎛ σ ⎞ 9 ⎛ σ ⎞6 ⎛ σ ⎞ ⎪ ⎪ 4⎢⎜ ii ⎟ − ⎜ ii ⎟ − ⎜ ii ⎟ + ⎜ ii ⎟ ⎥ r < rc ⎝r ⎠ f (r ) = ⎨ ⎢⎣⎝ r ⎠ ⎝ rc ⎠ ⎝ rc ⎠ ⎥⎦ ⎪ ⎪0 r ≥ rc ⎩
Taking a derivative of G with respect to εi+1 yields the following expression for updating the ε value at each iteration
(9)
If we used Ui(r) = εi f(r) in the ith iteration, the objective was to find εi+1 that provides the best fitting gi+1(r) for the g(r), such that the difference between Ui+1(r) = εi+1 f(r) and Ui(r) = εi f(r) approaches zero. In practice, we minimized the function G=
=
∫
2 ⎡ ⎛ g (r ) ⎞⎤ i ⎟⎟⎥ dr w(r )⎢Ui + 1(r ) − Ui(r ) − αkBT ln⎜⎜ ⎢⎣ ⎝ g (r ) ⎠⎥⎦
∫
2 ⎡ ⎛ g (r ) ⎞⎤ i ⎢ ⎥ ⎜ ⎟ w(r ) (ϵi + 1 − ϵi)f (r ) − αkBT ln⎜ ⎟ dr ⎢⎣ ⎝ g (r ) ⎠⎥⎦
(11)
εi + 1 = εi + αkBT
gi(r )
( ) dr
∫ w(r )f (r ) ln 2
g (r )
∫ w(r )f (r ) dr
(12)
Because the RDF is a discrete function, we used Simpson’s rule to handle the integral in eq 12. α is set to unity for the first iteration, and it is decreased by a factor of 0.8 after each iteration until the εi value is converged (with tolerance of 10−6).
■
RESULTS AND DISCUSSION Intermolecular Nonbonded Interaction Parameters. We parametrized our CG force field at two different temperatures: 25 and 50 °C. 50 °C is the typical gelation temperature for a dilute MC solution (
