Solvent-Driven Preferential Association of Lignin with Regions of

Aug 27, 2013 - The precipitation of lignin onto cellulose after pretreatment of lignocellulosic biomass is an obstacle to economically viable cellulos...
0 downloads 7 Views 5MB Size
Article pubs.acs.org/Biomac

Solvent-Driven Preferential Association of Lignin with Regions of Crystalline Cellulose in Molecular Dynamics Simulation Benjamin Lindner,*,†,‡ Loukas Petridis,*,† Roland Schulz,*,†,‡ and Jeremy C. Smith*,†,§ †

UT/ORNL Center for Molecular Biophysics, Oak Ridge National Laboratory, 1 Bethel Valley Road, Building 6011, Oak Ridge, Tennessee 37830-6309, United States ‡ Genome Science and Technology, University of TennesseeKnoxville, F337 Walters Life Sciences, 1414 Cumberland Avenue, Knoxville, Tennessee 37996, United States § Department of Biochemistry and Cellular and Molecular Biology, University of Tennessee, M407 Walters Life Sciences, 1414 Cumberland Avenue, Knoxville, Tennessee 37996, United States S Supporting Information *

ABSTRACT: The precipitation of lignin onto cellulose after pretreatment of lignocellulosic biomass is an obstacle to economically viable cellulosic ethanol production. Here, 750 ns nonequilibrium molecular dynamics simulations are reported of a system of lignin and cellulose in aqueous solution. Lignin is found to strongly associate with itself and the cellulose. However, noncrystalline regions of cellulose are observed to have a lower tendency to associate with lignin than crystalline regions, and this is found to arise from stronger hydration of the noncrystalline chains. The results suggest that the recalcitrance of crystalline cellulose to hydrolysis arises not only from the inaccessibility of inner fibers but also due to the promotion of lignin adhesion.



INTRODUCTION Lignocellulosic biomass is a highly complex material containing cellulose fibers laminated by other biopolymers, such as lignin and hemicellulose.1,2 The cellulose in higher plant cells is an abundant potential renewable feedstock for the production of biofuels, such as ethanol, via hydrolysis and fermentation.3,4 However, biomass is naturally recalcitrant to deconstruction.3 The industrial process for the production of biofuels normally involves pretreatment to disrupt cell-wall structure and increase yield.4 Pretreatment is followed by enzymatic hydrolysis which produces oligosaccharides as substrates for subsequent fermentation into ethanol.4 Recalcitrance itself is expressed in terms of the resistance to the hydrolytic activity of enzymes on cellulose and is usually determined by sugar yields. Several pretreatment techniques have been examined,5 each producing a distinct effect on the physicochemical features of biomass.6 One major thermochemical technique employs grinding, dilute acid, and high temperatures,7 resulting in pretreated biomass consisting largely of only cellulose and lignin. After dilute-acid pretreatment, lignocellulosic biomass still exhibits significant recalcitrance to hydrolysis.3,4,8 This has been attributed in part to the structure of the cellulose fibers, with amorphous forms being more hydrolyzable than crystalline.9 Furthermore, lignin remains in the sample and inhibits enzymatic activity. Although the molecular mechanism by which lignin has this inhibitory effect is not known in detail, there is evidence that it reprecipitates and binds both to the cellulose substrate10−12 and to the hydrolytic enzymes.13−15 A recent combination of small-angle neutron scattering (SANS) and molecular dynamics © 2013 American Chemical Society

(MD) simulation found model softwood lignin aggregates to have a rough, highly folded surface with pores capable of accommodating cellulolytic enzymes.16 Clearly, recalcitrance is strongly influenced by the physicochemical properties of lignocellulosic biomass.8 Among the available experimental techniques that have been used to probe biomass structure are X-ray scattering, electron12 and atomic17 force microscopy, SANS,18 NMR,19 and Raman spectroscopy.20 Considerable attention has been placed on cellulose, whose chains are crystalline to a varying degree of quality and quantity.21,22 The structures of two common naturally occurring crystalline forms of cellulose, I-β and I-α, have been characterized at atomic resolution experimentally.23,24 Other forms are categorized as paracrystalline21 or amorphous, denoting ultrastructures with distorted or disrupted interchain hydrogen bonding patterns, respectively. X-ray powder diffraction experiments provide a means of estimating bulk crystallinity indices,25 and microscope techniques have allowed local properties to be measured, such as cellulose fiber size distributions and pore sizes.11,26 However, the complexity of biomass has so far precluded detailed experimental characterization of the interaction of lignin with itself and cellulose at the molecular scale and, in particular, the effect of cellulose crystallinity on lignin reprecipitation is not known. Received: March 28, 2013 Revised: August 23, 2013 Published: August 27, 2013 3390

dx.doi.org/10.1021/bm400442n | Biomacromolecules 2013, 14, 3390−3398

Biomacromolecules



Atomic-detail information on molecular processes in model biomass systems can be provided by MD simulation. Previous MD work have studied the structure and dynamics of cellulosic oligomers,27 whole cellulose fibrils,28−31 and lignin.32 Early MD studies of cellulose−lignin association examined the binding of 10- and 20-unit guiacyl oligomers to cellulose,33,34 and these studies found that the adsorption of lignin onto different surfaces of a crystalline cellulose model yielded similar interaction energies. However, to study the aggregation processes of a statistically significant number of lignin molecules with a cellulose fiber, the required time and length scales are considerable and have only recently been rendered attainable for all-atom MD simulation, facilitated by the establishment of peta-scale supercomputing facilities and the formulation of specialized parallelization strategies.35 The requirements for a realistic simulation model for the present problem are discussed in Methods. The present study, which is one of the largest-ever atomistic biomolecular simulations to date, provides details on the processes involved in lignin precipitation onto cellulose and the dependence on cellulose crystallinity. Lignin shows a strong tendency toward aggregation with lignin and cellulose. The simulations reveal that the average direct interaction energy between lignin and cellulose is independent of cellulose crystallinity. However, lignin aggregation onto noncrystalline cellulose is less favorable than for the crystalline form, owing to the significantly different solvation properties of crystalline and amorphous cellulose. Hence, two factors strongly influencing biomass recalcitrance, cellulose crystallinity and lignin precipitation onto cellulose, are found to be interdependent.

Article

METHODS

Simulation and Modeling. Due to the heterogeneous nature of lignin the construction of a realistic model for lignin reprecipitation onto cellulose requires a large set of individual lignin molecules. Lignin clusters found experimentally range in size from 5 nm up to 10 μm,11,12 imposing a limit on the minimum size of the simulation box. Also, the degree of polymerization for cellulose is typically very large, ranging from 100 to 10000.36 The size of the system determines the time scale simulated, because lignin molecules must be able to translate within the simulation box, and the distance they travel is inherently diffusion limited. For a spherical particle with a radius of 2 nm, the self-diffusion constant provided by the Stokes−Einstein equation is on the order of 10−6 cm2 s−1, which results in an estimated translational diffusion of 17 nm in 500 ns. This determines the approximate minimum simulation time scale required to study lignin precipitation since lignin molecules must undergo multiple binding events with cellulose and other lignin molecules. Each lignin molecule modeled consists of 61 monomers, corresponding to a molecular weight of 13 kDa, that is, within the experimentally determined range.47 The monomer chemical composition and linkage were also obtained from experiment.48 The cellulose fiber model consists of 36 chains with a chain length of 160 monomers and is depicted in Figure S1A. The atomic starting coordinates for the cellulose are based on the Iβ crystal phase, derived from crystallographic studies,24 and the consensus model for the fiber.26 The noncrystalline cellulose model was generated by simulating crystalline cellulose at a high temperature (650 K) for 1 ns, which is well above the melting point for the internal hydrogen bonding network. The crystallinity of the cellulose was assessed by computing the 1D and 2D Bragg scattering diagrams49 and is shown in Figures S1B and S1C. The fully crystalline model exhibits strong Bragg peaks typical for Iβ, which are absent in the noncrystalline model. All models were hydrated with a 1 nm shell of explicit water, resulting in simulation sizes of 3.31, 3.43, and 3.80 million atoms for the NC, FC, and FN models, respectively (see Figure 1A and Results for model definitions). All models contain 1.56 MDa of biomass.

Figure 1. (A) Schematic of the initial configuration of the simulations. For each model, 52 lignin molecules were placed around the central fiber at distances randomly chosen within the interval of 0.2 to 1.0 nm and 0.4 to 2.0 nm for the NC model and the FC and FN models, respectively. The directional X-ray scattering diagram along the fiber axis is shown for both the crystalline and the noncrystalline cellulose, with the crystalline cellulose exhibiting strong Bragg peaks, which are absent in the noncrystalline model. (B) Visual representation of the initial and the final configurations of the systems. Lignin molecules belonging to the same cluster at the end of the simulation are colored identically. The lignin molecules participating in any given cluster tend to be in local proximity at the beginning of the simulation. 3391

dx.doi.org/10.1021/bm400442n | Biomacromolecules 2013, 14, 3390−3398

Biomacromolecules



The size of the cellulose fiber and the amount of lignin used in the simulation matches the experimentally observed ratio of molecular weights in softwood (cellulose/lignin = 1:0.52).48 The simulations were performed with the GROMACS MD suite50 using the TIP3P water model51 and the CHARMM Carbohydrate Solution Force Field (CSFF)52 for cellulose and the CHARMM Lignin Force Field.53 Since this work was completed, the carbohydrate force field has been revised.54 However, the changes made, which affect principally details of crystal structure geometries, are unlikely to impact any of the properties examined is the present article, which concern principally solvation and interfacial energies. Each of the models was simulated for 750 ns at 300 K using the NPT ensemble at 1 atm, with the NC model simulated twice with different initial starting velocities to examine the dependence of the aggregation on random fluctuations in the thermostat and barostat, giving a total of 3 μs simulation time. No significant difference was found between the two NC simulations (see Supporting Information), and hence, only data for one NC model are shown. The trajectories were saved for analysis every 5 ps. Further details on the molecular models and simulation protocol are provided in the Supporting Information. Aggregation Properties. Three different quantities were calculated to track lignin aggregation: the number of atomic contacts, N, the solvent accessible surface area, SASA, and the intermolecular interaction energy, E. The number of molecular contacts, M, and the molecular coordination numbers were derived from N. The above three quantities were found to be strongly correlated. The presence of the two different types of polymer, cellulose and lignin, requires the distinction between two different but coupled aggregation processes: lignin−lignin and cellulose−lignin. N, SASA, and E each form a 53 × 53 matrix at any given time step of the simulation, the first element representing the cellulose and the remaining elements corresponding to the 52 lignin molecules. This matrix is sparse since all measures are nonzero only for molecule pairs which are in spatial proximity. This approach enables each quantity to be pairwise decomposed. The calculation of N was performed with the tool g_mindist of the GROMACS MD suite.50 Two atoms were considered to form a contact if separated by 3 Å or less. The SASA was calculated with the VMD software55,56 employing a probe radius of 0.14 nm. A measure complementary to SASA is the buried surface area, BSA, which is the solute−solute interface not accessible to the solvent. The BSA of a solute molecule was calculated by subtracting its SASA from its total surface area, that is, BL(t) = S0L − SL(t) and BC(t) = S0C − SC(t), where S and B correspond to SASA and BSA, the suffixes L and C indicate lignin and cellulose, and S0 is the total SASA in the absence of L−L or C−L contacts. The large size of the trajectory turned the analysis into a high performance computing problem, requiring careful consideration of algorithmic design and data alignment, which is explained in the Supporting Information.

Article

RESULTS

The modeling approach is visually summarized in Figure 1A. Three models were simulated, each consisting of a cellulose fibril surrounded by 52 softwood lignin molecules in explicit water. Details of the molecular models are provided in Methods. The differences between the models lie in the initial random lignin placement and the crystallinity of the cellulose fiber and are visualized in Figure 1B. The following nomenclature is used: model NC model FC model FN

places the lignin molecules initially “near”, i.e., 2−10 Å, from a crystalline fiber places the lignin molecules initially “far”, i.e., 4−20 Å, from a crystalline fiber places the lignin molecules initially “far”, i.e., 4−20 Å, from a noncrystalline fiber

All lignocellulose simulation models showed the same overall behavior involving lignin molecules aggregating with each other and binding to the cellulose fiber. The precipitation of lignin onto cellulose and lignin self-aggregation can be quantified by the time-evolution of the total number of atomic contacts, N, for each of the three simulation models. N exhibited strong linear correlation with the buried surface area (Figure S5), and therefore, these quantities can be interchanged without affecting the conclusions. Cellulose−Lignin Aggregation. In all simulations, the number of cellulose−lignin contacts, N(CL), rises sharply over the first 100 ns before converging by the end of the simulations (Figure 2A), with a variation of less than 20% in the second half of the simulation. N(CL) is clearly different for each model throughout the whole simulation time. As will be shown later, this difference arises from the particular aggregation properties of the molecules (stickiness) and is not due to insufficient translational diffusion of the lignin molecules. Also, as may be visible from Figure 1, the edges do not play a significant role in the aggregation processes. Any edge effects should be considered an artifact of the limited size of the simulation, because experimental cellulose fibers are considerably longer. For the near-lignin/crystalline-cellulose (NC) model the near placement of the lignin molecules leads to a larger N(CL) at all times of the simulation than in the “far” models, in which the lignin molecules were initially placed further away from the fiber. Interestingly, a significantly higher number of cellulose-lignin contacts is seen for

Figure 2. (A) Total number of cellulose−lignin atomic contacts, N(CL). B) Total number of lignin−lignin atomic contacts, N(LL). For both graphs the inset shows the first 10 ns of each simulation. 3392

dx.doi.org/10.1021/bm400442n | Biomacromolecules 2013, 14, 3390−3398

Biomacromolecules

Article

far-lignin/crystalline-cellulose (FC) than far-lignin/noncrystalline-cellulose (FN), even though the lignin molecules were initially placed in a similar manner around the fiber. Thus, Figure 2A indicates that two factors affect the proportion of lignin reprecipitation onto cellulose: the average initial distance from the cellulose and the cellulose crystallinity. The slope of the number of atomic cellulose-lignin contacts is consistent with a decay equation, governed by two relaxation times (details provided in section 2.1 of the Supporting Information). The fast process has an associated relaxation time of about 10 ns for NC and 30 ns for FC and FN, which matches the observed rate at which new lignin molecules attach to the cellulose surface. The slower process can be interpreted as local reorganization of the lignin molecules on the cellulose surface with an estimated relaxation time of 150 ns for the NC and about 300 ns for FC and FN. The larger relaxation times for FC and FN indicate that an increase in initial cellulose−lignin distance results in a slow down of the reprecipitation process. Extrapolating the fitted functions to infinite time predicts that the amount of cellulose−lignin contacts converge against different values for crystalline and noncrystalline cellulose systems for time-scales significantly longer than simulated here. To obtain a thermodynamic understanding of the differences in lignin reprecipitation onto crystalline and noncrystalline cellulose, the cellulose−lignin and lignin−lignin interaction energy densities, ε, were calculated. ε is derived by dividing the interaction energy, E, by B, the buried surface area between cellulose−lignin or lignin−lignin. As is shown in Figure S5B,D, B and E exhibit strong linear correlation. No significant dependence of the cellulose−lignin interaction energy on cellulose crystallinity is found (Table 1). This is consistent with earlier

Table 2. Solvent−Biomass Interaction Constants solvent interaction energy densitiesa,b

ε

εcoul

εLJ

cellulose/lignin

NC FC FN NC FC FN

−47.65 ± 2.0 −51.80 ± 2.0 −50.33 ± 2.0 −52.66 ± 0.67 −52.95 ± 0.67 −52.02 ± 0.67

−20.43 −22.67 −24.10 −22.19 −22.24 −21.54

−27.22 −29.13 −26.24 −30.47 −30.72 −30.48

lignin/lignin

εcoul

εLJ

lignin crystalline cellulose noncrystalline cellulose

−74.1 ± 9.5 −94.0 ± 1.7 −107.3 ± 0.7

−63.0 ± 8.5 −83 ± 1.8 −102 ± 0.7

−11.2 ± 1.1 −10.0 ± 0.2 −5 ± 0.2

Energies are kJ/mol/nm2. bWater-solute constants are based on the simulation of isolated molecules. The standard deviation for water− lignin reflects the spread across all of the 52 lignin molecules.

noncrystalline cellulose. It was also found that, although the average number of hydrogen bonds made by cellulose (internal + solvent) per glucose molecule is very similar for the two models, with about 3.4 ± 0.1 for crystalline and 3.6 ± 0.1 for noncrystalline cellulose, the ratio of internal to solvent is markedly different, being 68 to 32% for the crystalline form compared to 47 to 53% for noncrystalline cellulose. The increased interaction energy between noncrystalline cellulose and water suggests that noncrystalline cellulose is effectively more hydrophilic than the crystalline form and explains why the FN model exhibits less lignin precipitation on to cellulose than the FC model (Figure 2A). This finding is further supported by the observation that lignin preferentially binds to the hydrophobic surfaces of crystalline cellulose, as shown in Figure 3, where the number of atomic cellulose−lignin contacts was decomposed into contacts with hydrophobic and hydrophilic cellulose surface. Based on the value for cellulose−lignin contacts for the hydrophilic chains of crystalline cellulose and all chains of noncrystalline cellulose are comparable, it is the hydrophobic surface of cellulose, which drives the increased reprecipitation of lignin onto crystalline cellulose. Another indicator for the hydrophobicity of cellulose was also investigated and is explained in detail in the “Hydrophobicity and Compressibility” section of the Supporting Information, where the compressibility was used as a dynamic probe for hydrophobicity.37−41 The lignin−water interaction energy densities are broadly distributed, ranging between −58 and −95 kJ/mol/nm2, and also were found to exhibit a strong correlation with the number of lignin−water hydrogen bonds but only a weak correlation with the total hydrophobic fraction of the surface area (Figure 4). An example of two lignin molecules with the same solvent accessible surface area and chemical topology but a very different number of hydrogen bonds with water and, thus, lignin−water interaction energy density is also given in Figure 4. In contrast to cellulose, the internal and solvent hydrogen bonds of lignin do not systematically compensate, and thus, lignin molecules exhibit a varying degree of unsatisfied hydrogen bonding groups. During the simulation, the individual lignin molecules undergo various conformational changes, which can be quantified by monitoring the difference in radius of gyration and is shown in Figure S3 in the Supporting Information, together with the difference in asphericity. The average radius of gyration tends to decrease slightly from 0 to 500 ns, which may be an effect of compactification during the aggregation. However, the large spread in radius of gyration changes renders this effect negligible, indicating that the average geometric properties of the lignin ensemble do not change during the reprecipitation. Lignin Self Aggregation. During the initial phase of lignin self-aggregation in the simulations, lignin molecules bind to

interaction energy densitiesa,b model

ε

a

Table 1. Lignocellulose Interfacial Interaction Constants interface

solute

a

Energies are kJ/mol/nm2. bThe assigned error bar is the uncertainty derived from the difference of the values between the two simulations of the NC model (second NC simulation: ε(CL) = −50.46 kJ/mol/nm2 and ε(LL) = −51.71 kJ/mol/nm2). The true error of εtotal depends on the variations in the possible aggregation pathways for each model.

findings suggesting that the lignin monomer binding interaction energy is independent of the molecular structure of cellulose.33,34 As a next step, solvation was investigated as a possible origin of the difference between the FC and FN models. To characterize the solute−water interactions, separate simulations were performed in explicit water of individual lignin molecules and of a crystalline and a noncrystalline cellulose fiber. The solute− water interaction energy was then calculated as the interaction energy per unit solvent accessible surface area ε = E/S. Table 2 shows that the interaction energy per unit area between cellulose and water is significantly larger for the noncrystalline than for the crystalline form. The difference in interaction energy density was traced back to differences in the capacity of the cellulose to hydrogen bond with water, which increases from 3.2 ± 0.1/nm2 for crystalline to 3.7 ± 0.1/nm2 for 3393

dx.doi.org/10.1021/bm400442n | Biomacromolecules 2013, 14, 3390−3398

Biomacromolecules

Article

Figure 3. Average number of atomic cellulose−lignin contacts per chain for three groups of cellulose surface chains: All, comprising all cellulose surface chains; Phobic, containing cellulose surface chains classified as hydrophobic; Philic, comprising surface cellulose chains classified as hydrophilic. The same assignment was used for the NC, FC, and FN models.

Figure 4. (A) Total SASA in relation to the fraction of hydrophobic for each lignin molecule. The linear regression yields a correlation coefficient of R2 = 0.714. This suggests that the fraction of hydrophobic atoms at the surface is size-dependent. (B) Lignin−water hydrogen bond density, which is the number of hydrogen bonds divided by the SASA, in relation to the interaction energy density for each lignin molecule. The strong correlation with a coefficient of R2 = 0.9872 suggests that the interaction energy density is governed by the hydrogen bonding capacity. (C) Molecular configuration for two lignin molecules with very different interaction energy densities. Both lignin molecules have approximately the same SASA. However, ID 34 has a significantly larger hydrogen bonding capacity with water. The color scheme we used is based on the element identity of each atom. Hydrogen atoms are colored in white, carbon is cyan, and oxygen is red.

molecular contacts. Molecular contacts are defined here as two molecules having at least one atomic contact, and the state of a lignin molecule is then defined by the number of molecular contacts formed. The states were labeled L for isolated lignin, LX for a lignin molecule in contact with X − 1 other lignin molecules and a leading C, if the lignin forms an additional contact with cellulose. The model and the time-dependent population for all states are provided in Figure 5. At 50 ns, FN shows a striking increase in the L3 and L4 states, indicating the formation of larger lignin clusters without cellulose contact, which coincides with a strongly reduced population of the CL2 state, which corresponds to two lignin molecules bound to each other and to cellulose. Hence, the reduced affinity of noncrystalline cellulose for lignin during the initial 50 ns leads to a long-term increase in the lignin self-aggregation.

macromolecules in their local vicinity. As is apparent from Figure 2B, during the first 50 ns all models show a similar amount of lignin−lignin aggregation, which is a direct consequence of the models having similar initial distances between the lignin molecules. However, after 50 ns, the FN model diverges from NC and FC, exhibiting a higher degree of lignin self-aggregation. Hence, noncrystalline cellulose also promotes lignin selfaggregation. As in the case for cellulose−lignin, the slope of the number of atomic for lignin−lignin contacts in Figure 2B is consistent with a decay equation. However, only one relaxation time was required to achieve an adequate fit, with NC at about 140 ns, and FC and FN at about 180 ns. To further quantify the self-aggregation effect, the lignin aggregation process was decomposed into a discrete state model, in which each state represents a distinct arrangement of 3394

dx.doi.org/10.1021/bm400442n | Biomacromolecules 2013, 14, 3390−3398

Biomacromolecules

Article

persistent binding of lignin to cellulose or another lignin molecule. We define the stickiness, S(t), as the conditional probability of finding a lignin molecule bound to another molecule at time t, given that the two molecules formed a contact at time ≤ t (formal definitions are provided by eqs S3 and S6). Cellulose−lignin and lignin−lignin stickiness, SCL(t) and SLL(t), respectively, are shown in Figure 6A,B. SLL(t) converges to 0.70 ± 0.05 for all models within the first 200 ns, which means that on average each lignin molecule retains 70% of the molecular contacts it initiates with other lignins. The association of a given lignin molecule with cellulose is influenced by the intrinsic propensity of the two molecules to bind in aqueous solution and also by extrinsic factors, such as the presence of other lignins that may block access to the cellulose. To probe the intrinsic factors contributing to stickiness of lignin on cellulose, SCL(t) was computed for the subset of lignin molecules which enter the CL state first, that is, bind first to cellulose before binding to other lignins, which makes the derived SCL(t) less sensitive to competing binding of lignin to cellulose. As shown in Figure 6A, SCL(t) converges to 0.80 ± 0.01, 0.75 ± 0.05, and 0.67 ± 0.01 for NC, FC, and FN, respectively. The difference in SCL(t) between the crystalline and noncrystalline cellulose models is consistent with the overall reduced lignin−noncrystalline cellulose association. Furthermore, the large values of SCL(t) are quantitatively indicative of the high degree of lignin stickiness for cellulose and explains why the initial conditions in the simulation have a significant effect on the morphology of the final metastable states. The large value of the lignin−lignin stickiness at 750 ns leads to the conclusion that each system has converged to a metastable state, because it shows that the lignin molecules have not performed a significant number of association and disassociation events. The latter is true because the stickiness is a time-dependent quantity and formally converges against zero for infinite times as lignin molecules come eventually in contact with all other lignin molecules. Here, in the limit of infinite simulation time, SLL(t) is expected to converge against the range between 0 and 0.1 for an equilibrated system, with the maximum value determined by the highest aggregation state for lignin (5), divided by the number of lignin molecules in the simulation (52). The differences between the models NC and FC are an indicator for the glassiness of the system, that is, the inability of lignin bound to cellulose to undergo large-scale molecular

Figure 5. Time dependence of the populations of the molecular lignin coordination number shown for each model and state.

Metastability and Stickiness. Because NC and FC were simulated with the same structural models for cellulose and lignin, the quantitative difference in N(CL) between NC and FC provided in Figure 2A indicates that the simulations have reached different metastable states. One reason for this metastability is the “stickiness” of lignin for both cellulose and lignin, that is, the

Figure 6. (A) Time dependence of the stickiness for lignin−lignin and (B) cellulose−lignin molecular contacts. Cellulose−lignin stickiness shown only contains lignin molecules which enter the CL state first, i.e., bind to cellulose before to any other lignin, which emphasizes the intrinsic binding propensity of lignin with cellulose. 3395

dx.doi.org/10.1021/bm400442n | Biomacromolecules 2013, 14, 3390−3398

Biomacromolecules

Article

Figure 7. (A) Probability distributions of atomic contacts between cellulose and lignin based on cluster morphology, averaged over 50k frames (250− 500 ns, minor structural rearrangements with a variation in total atomic contacts of less than 20%). (B) Bubble iagram of the Probability (Size of Bubbles) of lignin clusters of a given size (x axis) and their cellulose association (y axis), quantified by the fraction of lignin molecules in contact with cellulose. Cluster sizes of 1 mark lignin molecules in the L or CL state. FN has large clusters with low association, while NC results in medium sized clusters with high cellulose association. For visualization, the data for NC was shifted by 0.25 on the x axis to the lower values and FC was shifted 0.25 to higher values. (C) Probability distributions of atomic contacts between lignin molecules with surrounding lignin, averaged over 50k frames (250−500 ns).

Figure 8. Pairwise cellulose−lignin interaction energy as a function of buried surface area. Each point represents a group of lignin molecules in the same aggregation state (CL, CL2, CL3, CL4, CL5). The underlying distribution is provided in Figure S6 of the Supporting Information. The data was calculated from 2000 configurations between 740 and 750 ns of simulation time.

number of contacts for FN (280) is less than the maximum for FC (400). Similar information can be retrieved from Figure 7B, which shows the probability of a lignin cluster plotted against the cluster size and the degree of association with cellulose (the number of molecular lignin−cellulose contacts) and reveals, for example, that the FN model contains a large lignin cluster (17 molecules) with hardly any cellulose association (