New Mechanistic Insights into the Lignin β-O-4 Linkage Acidolysis with Ethylene Glycol Stabilization Aided by Multilevel Computational Chemistry.

Acidolysis in conjunction with stabilization of reactive intermediates has emerged as one of the most powerful methods of lignin depolymerization that leads to high aromatic monomer yields. In particular, stabilization of reactive aldehydes using ethylene glycol results in the selective formation of the corresponding cyclic acetals (1,3-dioxolane derivatives) from model compounds, lignin, and even from softwood lignocellulose. Given the high practical utility of this method for future biorefineries, a deeper understanding of the method is desired. Here, we aim to elucidate key mechanistic questions utilizing a combination of experimental and multilevel computational approaches. The multiscale computational protocol used, based on ReaxFF molecular dynamics, represents a realistic scenario, where a typical experimental setup can be reproduced confidently given the explicit molecules of the solute, catalyst, and reagent. The nudged elastic band (NEB) approach allowed us to characterize the key intermolecular interactions involved in the reaction paths leading to crucial intermediates and products. The high level of detail obtained clearly revealed for the first time the unique role of sulfuric acid as a proton donor and acceptor in lignin β-O-4 acidolysis as well as the reaction pathways for ethylene glycol stabilization, and the difference in reactivity between compounds with different methoxy substituents.

Introduction

The highly abundant aromatic biopolymer lignin is a key renewable resource for the production of important aromatic chemicals,[1−4] which led to the development of proficient depolymerization strategies to obtain valuable monophenolic products.[5−9] Among these, acidolysis(10−13)—one of the most classical methods for breaking down the structure of lignin—has re-emerged as a viable way of lignin depolymerization, especially in combination with robust stabilization strategies,[14−16] which can suppress undesired recondensation phenomena and markedly increase monomer yield.[17]

Elegant studies have established the main principles of the scission of the β-O-4 moiety in lignin and model compounds,[10,12,13,15,16,18−28] and recent computational work allowed detailed mechanistic insight to be gained through these processes.[29−32] It has been generally recognized that the scission of the β-O-4 moiety starts with protonation and dehydration of the α-OH, readily furnishing a benzylic carbocation, which subsequently undergoes two productive C–O cleavage events (Figure A) (a) yielding the family of Hibbert ketones via rehydration (C3-pathway) or (b) leading to the formation of C2-aldehydes via deformylation/rehydration (C2-pathway), whereby the extent of the C2 vs C3 pathways depends strongly on the reaction conditions and the type of mineral acid used.[11,12,13,15,20,33−38] It has to be noted that the benzylic carbocation also engages in nonproductive condensation phenomena leading to recalcitrant C–C bonds (Figure A). Over the past few years, several computational investigations focused on the elucidation of the mechanisms of lignin depolymerization in homogeneous[35−37,39−49] or heterogeneous phases.[50−55] The general approach is to consider gas-phase structures of reactants and products, eventually immersed in a continuum dielectric medium, and to identify transition states from several possible pathways suggested by chemical intuition and experiments. Among these excellent investigations, in a recent experimental–computational work, Beckham and co-workers[29,30] employed state-of-the-art density functional theory (DFT) calculations with an appropriate choice of the functional, a well-balanced combination of explicit solvent molecules with a continuum solvent description (tuned cluster/continuum methods) and microkinetic modeling to investigate in detail the nature of these C–O scission events. They have elaborated on various possible acidolysis mechanisms (without stabilization) and identified the importance of using realistic aromatic substitution patterns in β-O-4 model compounds, as this greatly influences the rate of C–O bond scission and condensation. Significant rate differences could be assigned to the formation of key intermediates leading to the aldehyde and ketone cleavage products. Their success in finding satisfactory agreement between experimental and modeling was due to an exhaustive exploration of the complex structures, strategical incorporation of the quantum chemistry results with a microkinetic model, and an optimal calibration of the parameters.

Figure 1

(A) Schematic representation of lignin β-O-4 cleavage in lignin and major pathways involved. (B) Acidolysis reaction scheme for the β-O-4 cleavage in lignin with ethylene glycol stabilization. (C) Approach adopted in this work: A combined experimental and multilevel computational strategy

By simple and advanced model compounds[15,16] and lignin studies,[56,57] we have introduced the importance of stabilization of unstable C2-aldehydes formed upon acidolysis by trapping these using diols to obtain more stable C2-cyclic acetals (Figure B), as a means to prevent condensation phenomena and markedly increase useful monomer yield. We found that in aprotic solvents, catalytic amounts of acid were sufficient for a rapid β-O-4 cleavage reaction, while ethylene glycol served particularly efficiently in the formation of 1,3-dioxolanes, and demonstrated that under certain reaction conditions (EG and acid concentration), incorporation of EG into the benzylic position took place, which also stabilized the previously mentioned benzyl carbocation intermediate (Figure B).[58] Very recently, we presented a novel, diol-assisted-fractionation protocol (DAF) of pine lignocellulose as a metal-free “lignin-first” strategy, applying EG, H2SO4, and the green solvent dimethyl carbonate (DMC) (Figure B).[58]

Motivated by these studies, here we wish to exploit the great potential offered by atomistic multiscale approaches in understaning fundamental questions related to the acidolysis of the β-O-4 moiety, and to gain new mechanistic insights into acidolysis in conjunction with our unique stabilization strategy. Specifically, we wish to gain mechanistic understanding regarding existing reaction pathways and intermediates—with particular attention to the use of EG as a stabilization agent, the possible influence of the type of acid and its counterion, anion–charge interactions, and solvent effects. To address all these complex and interwoven points, a combined experimental and multilevel theoretical approach is envisioned. In relation to the previous modeling studies, our multiscale computational protocol, based on ReaxFF Molecular Dynamics, represents a realistic scenario, where explicit molecules of the solute (model compounds 1-G/S, Figure C), catalyst, and stabilizer move freely and can interact/react with each other in the crowded environment that reproduces a typical experimental setup confidently. By reducing the dimension of the system from thousands of atoms to less than a hundred, the key processes identified at the ReaxFF level are then scrutinized and refined via quantum chemistry simulations. Reaction paths are characterized in detail through the NEB approach, disclosing all the reaction steps leading to key reaction intermediates and products (Figure C, 3-G/S, 4-G/S, 5-G/S, 6-G/S, 7-G/S)

The methodology developed here is capable of scrutinizing the role of H2SO4, DMC, and EG (stabilization agent) and shines a light on various proton transfer events involved in the cleavage of β-O-4 models. Thus, it provides significant insight into lignin acidolysis, which was so far not disclosed at the molecular level. The developed mechanistic insight will serve to improve lignin and lignocellulose processing and maximize the yield of biobased aromatics.

Results and Discussion

To elucidate the role of the various species (H2SO4/EG/DMC) during lignin acidolysis, we chose a combined experimental and computational approach for studying the reactivity of model compounds 1-G and 1-S that are simple representations of the internal β-O-4 linkage in lignin (Figure e)[26] flanked by G- and S- units and thus a good testbed for our modeling approach. Additionally, their use allowed us to target the C2-pathway directly without the formation of Hibbert ketones by the C3-pathway. Model compounds 1-G and 1-S were subjected to acidolysis conditions at 140 °C, in the presence of 5 mol % sulfuric acid (to 1-G or 1-S) in DMC with or without EG (Figure e, path 1 and 2) and the corresponding product formation profiles are shown in Figure . In the absence of EG, the formation of guaiacol or syringol 5-G/S and C2-aldehydes 6-G/S (Figure e, Path 1) were detected as 1-G/S acidolysis products, while 5-G/S and the corresponding C2-acetals 7-G/S (Figure e, Path 2) were formed in the presence of EG. In both cases, vinyl-ether 4-G/S was observed as an intermediate.

Figure 2

Acidolysis reaction profiles of (a) 1-G no EG (path 1); (b) 1-G with EG (path 2); (c) 1-S no EG (path 1); (d) 1-S in path 2 (with EG) as determined by HPLC/GC(FID). Reaction conditions: substrate (1-G or 1-S, 1 equiv) 0.016 mmol/mL in dimethyl carbonate (DMC), internal standard (1,2,4,5-tetramethylbenzene, 1.7 equiv), H2SO4 (5 mol % to 1-G or 1-S), 140 °C and for (b) and (d) ethylene glycol (EG, 4 equiv), nitrogen atmosphere.

During the acidolysis of 1-G, full conversion was reached in 15 min in the case of path 1 (no EG, Figure a) and in 1 h for path 2 (with EG, Figure b) with the initial formation of the vinyl ether 4-G in up to 50–60% selectivity, followed by the formation of C2 aldehyde 6-G was detected in 56% yield together with guaiacol 5-G (90%) after 2 h for path 1. On the other hand, with EG stabilization (path 2), C2-acetal 7-G and guaiacol 5-G were observed in the same, excellent 94% yield after 2 h as a result of efficient stabilization, in accordance with previously reported data.[15] Product formation profiles were measured at lower acid concentration (2 mol % instead of 5 mol %, see also Supporting Information (SI) Figure S6a and b), thus lower reaction rates as reflected by slower 1-G conversion, which allowed for better assessment of the formation and gradual disappearance of early reaction intermediates 3-G and 4-G that are formed prior to the actual β-O-4 cleavage reaction. In the absence of EG stabilization, only vinyl ether 4-G was observed as intermediate (Figure S6a), while with EG stabilization (Figure S6b), in addition to 4-G, the formation of α-EG benzyl ether 3-G was observed in the first 1.5 h of reaction up to around 17% after being consumed again. Interestingly, under these milder conditions, 6-G could be detected in small quantities (up to 2% after 2 h, Figure S6b), indicating either a very fast reaction with EG or the possibility of a different mechanism that may involve a direct pathway to 7-G without forming 6-G (vide infra).

The rate of conversion for the syringol based model 1-S turned out to be lower than 1-G both with and without EG (Figure d and c, respectively), reaching full conversion in 1 h (Figure c) and 2 h (Figure d), respectively. In path 1, the amount of aldehyde 6-S reached low yield even after 4 h (30% yield, Figure c), together with 55% 5-S. In path 2 (Figure d), comparable amounts of 5-S and 7-S were found (51% and 50%, respectively), indicating the necessity of efficient stabilization again. In the case of this more oxygenated model compound 1-S (Figures c and 2d), even though full conversion was reached after 4 h, it was possible to identify only a total of 60% of products (4-S 6%, 5-S 55%, and 6-S 30% for path 1; and 4-S 12%, 5-S 51%, and 7-S 50% for path 2), suggesting the possibility for condensation pathways that would result in higher Mw products or other nonvolatiles. A much better mass balance could be achieved by applying lower acid concentration (2 mol % instead of 5 mol %, Figure S6c and d) where 96% of the products could be quantified in the presence of EG (Figure S6d, 4-S 18%, 5-S 66%, 7-S 68%, 3-S 12% after 96h). In addition, under these conditions the intermediate 3-S (α-EG ether, probably not stable enough at 5 mol % acid and the selected EG concentration) and traces of 6-S were detected (Figure S6d). In path 1 (Figure S6c), 6-S was detected in 22% yield after 24 h which then disappeared to likely form condensation products. Furthermore, it should be noted that, especially in the low acid concentration regime (2 mol % sulfuric acid), 1-S conversion was found to be much slower than that of 1-G (up to 24 h to reach full conversion for 1-S while 45 min for 1-G). It is important to mention that in our system, DMC can also partly react. The formation of traces of ethylene carbonate (1,3-dioxolan-2-one) could be detected, but no carbonate modified 1-G/S or derived intermediates and products were found, as we also observed in a previously study of this reaction system in softwood.[58]

In summary, these results point to similarities but also significant differences between reactivity patterns of G and S based models with and without EG stabilization. The latter led to the formation of the more stable 7-G/S instead of 6-G/S, but also with 4-G/S as intermediate. Low acid concentration (2 mol %) allowed detection of 3-G/S in the presence of EG, and the low amount of 6-G/S formed under these conditions led us to consider an alternative route for the formation of 7-G/S (see below). Furthermore, there were clear differences in the conversion rates of 1-G and 1-S also in the presence of EG and at different acid concentrations. In addition to the catalytic H+ equivalents that H2SO4 offers, we also suspected a more direct role of the acid in the reaction mechanism, as previous literature clearly points at counterion effects in this reactions.[15,38,59] For example, H2SO4 is known to facilitate proton transfer reaction[60−63] of the type involved in the steps by which β-O-4 cleavage reaction can be explained.[32]

Thus, we aimed to further investigate the fundamentals of these observations involving reactions with H2SO4/DMC versus H2SO4/DMC/EG at the molecular level using multiscale modeling to elucidate intermolecular interactions among the various species and reaction mechanisms. More specifically, the acidolysis mechanism was examined comprehensively for paths 1 and 2, focusing on key catalyst–substrate interactions and underpinning specific molecular mechanisms that cause the different reactivity of 1-G and 1-S.

To create a picture of the environment surrounding the solute (1-G) under the applied acidolysis conditions, we started with the comparative analysis of the sampled configurations during the MD simulations with and without EG (Figure ). Initially, the solvation of the 1-G molecules was investigated, good solvation (Figure S1), but still, sufficient interaction between 1-G molecules forming aggregates of stacked and T-shaped structures (Figures and S2). These small clusters of two or three molecules were characterized by parallel or perpendicular orientations of the rings with compact or chain-like morphology. These types of interactions are essentially those that lead to expected aggregation in the more complex lignin structures and need to be overcome by a careful solvent selection.

Figure 3

(a) Percentage probability of finding X number of H2SO4 and EG molecules within 3.5 Å of the hydroxyl oxygen of 1-G in both simulations with and without EG. H2SO4 showed the same probability in both simulations. (b) Atom–atom radial distribution functions between the hydroxyl oxygen of 1-G and the carbonyl oxygen of DMC (red and blue lines), the oxygens of EG (green line), and the hydroxyl oxygens of H2SO4 (magenta and cyan lines).

Figure 4

Example of two 1-G molecules with stacked rings. Snapshot extracted from the ReaxFF MD simulations in a box containing H2SO4/DMC/EG molecules.

Looking further at the 1-G interactions, it was shown that the oxygen of the α-OH of 1-G (Figures and 2e) was frequently in contact with an H2SO4 molecule (80%) at about 2.7 Å (Figure a). This interaction was stronger than that with EG, which was found less frequently (50%) at close distance (2.8 Å) (due to its weaker interactions). Examination of the RDF profiles showed that EG competes with both the catalyst and the solvent, as evidenced by the broadening and shift of the H2SO4 peak to a longer distance (2.8 Å). The RDF trends also suggest that in the mixed solution, DMC has a slightly higher probability than EG of occupying the catalytic site and that the interaction of the α-OH with the sulfuric acid is dominant (Figure b).

Given these premises, we could speculate, in line with the literature, that the initial step of the reaction mechanisms is indeed the protonation of the α-OH by H2SO4 even in the presence of EG, followed by the release of a water molecule and the formation of carbocation 2-G (Figures and 2e), en route to 4-G.

Possible reaction paths, obtained with the NEB approach, to simulate the first step, intermediates, and final products (5-G and 6-G) without EG are shown in Figure . In the first step (Figure , step 1) 1-G undergoes an E1 elimination where the acid protonates the most basic atom in the molecule, a weakly basic water molecule is expelled, and a carbocation is formed (saddle point in the NEB profile). Then, the carbocation loses a proton in the β-position, which is released to the same formed hydrogen sulfate molecule to maintain neutrality and leads to the formation of 4-G. The cleavage of the C–H bond at the β-position is an intermediate minor reaction, whereas the formation of the carbocation intermediate is the rate-limiting step. This depends on the activation energy for breaking the connection between the carbon and the water molecule that is approximately 25 kcal/mol. This value perfectly agrees with the one previously observed for ΔG, for similar benzylic carbocation intermediates.[29,31]

Figure 5

Minimum energy path and molecular structures of the reaction mechanisms. Step 1: formation of 4-G starting from 1-G. Step 2a: formation of intermediate 8-G, where the hydroxyl moiety is connected to Cβ, starting from 4-G. Step 2b: formation of 6-G and 5-G starting from 8-G.

In the second step, the produced 4-G undergoes a subsequent H2SO4 mediated addition of a hydroxyl to the β position (8-G, Figure , step 2a). Interestingly, in our system, H2SO4 is shown to play the key role of proton-transfer bridge in all steps. According to the NEB profile (Figure , step 2a), a water molecule reaches the double-bond adjacent to Cβ; the acid donates a proton to Cα, and to maintain neutrality, the water releases a proton to the negatively charged hydrogen sulfate molecule resulting in the formation of the Cβ-OH. The transition structure is again a carbocation, but the energy barrier, which is around 15 kcal/mol, is much lower than the one estimated for the first step. The cleavage of the β-O-4 bond and the formation of 6-G and 5-G (Figure , step 2b) are activated by an acid proton interaction and subsequent transfer to the ether oxygen in the β-O-4 bond and a subsequent reprotonation of the acid by the −OH moiety (formally a protonated aldehyde at this stage). This mechanism is characterized by a shallow energy barrier of approximately 3.5 kcal/mol and a molecular configuration where the two transferred hydrogens are both connected to the produced fragments at an average distance of about 1.05 Å to the oxygen atoms interacting with two different oxygens within the sulfuric acid molecule.

The picture emerging from the different reaction mechanisms shown in Figure suggests that the route to the 5-G and 6-G is a multistep process occurring through the formation of a first carbocation intermediate (2-G) that is the rate-determining step due to its higher energy barrier. This intermediate is a high-polarity transition state and, as such, should be stabilized to improve its reactivity. A possible stabilization could be obtained with moderate/high-polarity solvents and also by including neutral/weak nucleophiles, as exemplified in a computational investigation focused on the acid-catalyzed hydrolysis of the lignin β-O-4 bonds in model compounds in an ionic liquid environment using DFT calculations and continuum models.[31]

These results obtained with ReaxFF MD simulations showed that H2SO4 has a crucial role in all steps functioning as proton donor and acceptor. To study this experimentally, we performed 1-G and 1-S acidolysis in previously applied conditions now using HCl (2 mol % to the substrate) as a catalyst (Table S1), and interestingly, no starting material conversion was observed. This is consistent with our previous results[15] where the lignin β-O-4 model compound 2-(2-methoxyphenoxy)-1-phenylethan-1-ol gave 0% conversion in toluene and 1,4-dioxane with HCl (10 mol %, 2 h). However, when we raised the HCl concentration (20 mol %, Table S1) we could observe the formation of 4-G/S in both cases (4% for 1-G and 10% for 1-S with 6% and 14% conversion, respectively). We speculate that this can be due to the presence of more than one molecule of HCl involved in the mechanism of donating and accepting protons at higher concentrations, whereas only one molecule is necessary in the case of H2SO4. This underscores the unique nature of H2SO4 action by an interplay of its strong acidity (and so its ability to perform the acidolysis reaction) and ionic size with multiple sites that can be involved in the donation and reception of protons, as discussed above. Considering the rate-determining step of 4-G formation, 4-G is formed after abstraction of the proton in the β-position of the carbocation intermediate by negatively charged hydrogensulfate. The latter step was widely discussed in the earlier acidolysis literature. Yokoyama and Matsumoto[24] reported the acidolysis of 1-G using HBr in 82% aqueous dioxane at 85 °C, and they proposed the β-proton abstraction by the solvent (water or dioxane). Nonetheless, in a later study using the C3 analogue model compound of 1-G, the participation of bromide anion in the general acidolysis mechanism was suggested, but not confirmed with calculation.[27] Overall, our results highlight that the interplay between the polarity of the solvent with its relative interaction with the catalyst, combined with the proton transfer capabilities of the catalyst, are important features for the lignin β-O-4 scission reaction to proceed as also demonstrated by Ekerdt et al.[35] or Matzumoto et al.[25] studying the model compound cleavage in acidic imidazolium-based ionic liquids or in 1,4-dioxane/water with HBr, H2SO4, and HCl, respectively. In our case, the cleavage reaction is particularly efficient in the presence of H2SO4 due to beneficial proton transfer interactions that facilitate key transition states that previous computational methodologies were not able to reveal.

To disclose the effects due to the combination of the molecular species considered above in the presence of EG (path 2), we reinvestigated in the reaction steps shown above, focusing on (a) the formation of 3-G (α-EG-incorporation) and (b) the formation of 7-G via different pathways and the role of EG and sulfuric acid herein.

The ethylene glycol adduct 3-G can be formed from 1-G via 2-G (Figure , equil. 1) or directly via a SN2 substitution mechanism. The obtained NEB profile of the formation of 3-G shown in Figure shows that this transformation goes via a similar carbocation intermediate (2-G) as does the formation of 4-G from 1-G (Figure , step 1) as the most common pathway speculated in the literature for substitution at the α-position of the β-O-4 with alcohols.[60−64] The energy barriers found of approximately 24 kcal/mol for the 1-G to 3-G and the 1-G to 4-G routes are also similar. In addition, the conversion of 4-G to 3-G was considered (Figure S3), and a 10 kcal/mol barrier was found showing that this reaction could also contribute to the formation of 3-G, which is typically not considered.

Figure 6

Minimum energy path and molecular structures of the reaction mechanisms. Equil. 1: the formation of 3-G starting from 1-G. Step 3: first, the formation of 9-G, where EG is connected to Cβ, starting from 4-G, followed by the formation of 7-G and 5-G directly from 9-G.

The formation of 7-G from 4-G can follow several pathways. The first considered option is that 6-G is formed according to the same mechanism described for pathway 1, which then undergoes acid-catalyzed acetalization with EG to give 7-G (Figure S4). This mechanism consists of two stages according to the general mechanism of cyclic acetal formation, namely, the formation of a hemiacetal (10-G, Figure S4a) first, and then the formation of the cyclic acetal by the expulsion of a water molecule (Figure S2b). The two steps are characterized by energy barriers of about 2.5 and 25 kcal/mol, respectively; the second step is much more difficult, as expected.

As an attractive alternative pathway, we considered the direct conversion of 4-G to 7-G and 5-G with the direct involvement of EG (Figure , step 3). Here, EG takes over the role of water as a nucleophile that adds to Cβ as shown in step 2 in Figure . This results in 9-G with an energy barrier of about 18 kcal/mol, which agrees with the energy range of 14–16 kcal/mol estimated to attach a hydroxyl moiety to those sites. Once EG has incorporated into the Cβ position, while the scission of the phenyl ether is activated by the protonation of the C–O bond and the simultaneous attack of the −OH moiety of EG accompanied by simultaneous proton transfer to the hydrogen sulfate anion. Similar to the formation of 6-G via step 1, the creation of a carbocation (NEB saddle point) is the rate-limiting step with an energy barrier of around 23 kcal/mol.

These results suggest that both mechanisms for the formation of 7-G are feasible (Figures ,6, and S2), but the direct alternative pathway has slightly lower barriers and is likely favored by the excess EG present in the reaction. This is supported by the experimental data, since 6-G was not observed when 5 mol % sulfuric acid was used (Figure b) and found in very low concentrations with 2 mol % acid (maximum 2% yield) when EG was present during acidolysis of 1-G (Figure S6b). According to the energy profiles determined for the pathway via the aldehyde formation, the largest barrier is found in the step going from 6-G to 7-G, indicating that 6-G would build up in solution, should this pathway be dominant.

Overall, the primary rate-determining step for both paths is the very first step (formation of 4-G from 1-G via 2-G), which is only mildly influenced by the competing pathway for the formation of 3-G in the presence of ethylene glycol. Thus, to explain the rate difference between 1-S and 1-G, especially observed at lower acid concentrations, we have focused on the first reaction step of forming 4-S via 2-S. Inspection of the NEB (Figure ) for 1-S to 4-S reveals that the barrier of the minimum energy path for step 1 (formation of the carbocation 2-S) is around 35 kcal/mol, which is significantly different compared to the corresponding energy barrier of 1-G to 4-G via carbocation 2-G (25 kcal/mol). A possible explanation was obtained by resorting to the partial atomic charges on the α-C–O (OH) bond extracted from the Natural Bond Orbital (NBO) analysis,[61] as suggested by Pelzer and co-workers.[30,65] Following their suggestions, we calculated the charge transfer at the saddle points by comparing the partial charge difference on the α-C–O in the reactants 1-G (1-S) and in the respective saddle structures, with the total charge transfer between reactants and products, after the appropriate normalization. We found that for 1-G, only 23% of the total charge that should be transferred during the reaction occurred before the saddle point, whereas in the 1-S case, the corresponding charge percentage transferred before the saddle point was 68%. This means that less charge needs to be transferred to obtain the transition state for 1-G than for 1-S, which clarifies the higher activation energy.

Figure 7

Minimum energy path and molecular structures of the reaction mechanism depicting the formation of 4-S starting from 1-S. In the inset, superimposed structures collected along the path.

Furthermore, a careful examination of the molecular configurations sampled along the path (inset in Figure ) suggests that major conformational rearrangements involving the acid and the methoxyl group of the ring nearby (not present in 1-G, Figure S5) took place. This also indicates that the action of the acid is sterically hindered by the added substituents and that the acid cannot find an optimal position to donate/abstract the hydrogens from the reaction location. These effects seem to be significant enough to counteract the enhanced inductive stabilization by the increase in the electron density of the aromatic ring due to the additional donating methoxy group present 1-S when compared to 1-G, thus leading to a significantly lower overall reaction rate. This is in line with what has been reported by Matsumoto, Yokoyama, and co-workers[66] which studied the formation rate of benzyl cation intermediate from p-hydroxyphenyl, guaiacyl, and syringyl nucleus using different model compounds. They found that the formation rate of syringyl nucleus benzyl cation intermediate decreased compared to the guaiacyl nucleus. In fact, the electron-withdrawing effect of the extra meta-substituting methoxy group present in the syringyl nucleus would destabilize the benzyl cation intermediate, decreasing its formation rate.

Earlier computational work using geometry optimization by energy minization on similar model compounds concluded that −OMe groups can play a role in stabilization of the carbocation by an intramolecular lone-pair interaction via specific structural arrangements.[35] Under our experimental conditions, however, this contribution seems to be less pronounced. This may be due to the specific solvent environment and/or considerations related to the subsequent β-H abstraction step. Our approach offers insight into the intermolecular interactions and the need for a favorable substrate orientation for an efficient β-H abstraction step, in which the steric hindrance of the −OMe group can explain the overall negative effect of this substituent on the reaction rate together with the observed charge transfer requirement.

Conclusions

A combination of a multiscale (from thousands of atoms to less than a hundred), multilevel (from classical simulations based on force fields to quantum chemistry calculations) theoretical description of all reaction steps responsible for the acidolysis of 1-G and 1-S to the final products (6-G/6-S or 7-G/7-S and 5-G/5-S) with experimental data was performed to elucidate the role of the different species present in the system (5 mol % H2SO4 in DMC) in the absence and presence of EG at 140 °C. Overall, the RDF trends showed that in the mixed solution of H2SO4/DMC the interaction of the α-hydroxyl oxygen of 1-G with the sulfuric acid is dominant in both cases (with or without EG) and the likely starting point of the reaction. The unique role of H2SO4 as proton acceptor and donor in this system was clarified and supported by experimental data. The involvement of EG in stabilization was computationally elucidated, suggesting two parallel pathways resulting in acetal 7-G. Here, the direct EG-addition to the vinyl ether intermediate 4-G in β-position to yield 9-G followed by direct cleavage to acetal 7-G seems to be preferred over the pathway via aldehyde 6-G that undergoes subsequent acetalization with EG. Additionally, the different reactivity of 1-G and 1-S, with an additional −OMe group on each aromatic ring, was explained in terms of the relative amount of charge transfer before the saddle points that in 1-S was three times higher than in the other case and steric impediments. The satisfactory results of this combined experimental and computational investigation demonstrate the power and efficiency of the work strategy and represent a step forward toward the definition of effective methods for rational reaction design. The object of future work will be to investigate the role of sulfuric acid in promoting the C2 versus C3 pathway which is fundamental to understand the behavior of lignin from different wood sources.

Materials and Methods

General Experimental Procedure for Lignin β-O-4 Model Compounds Acidolysis Reaction

Acidolysis reactions were carried out using previously reported conditions.[58] Two different procedures were followed: in one case, 1-G and 1-S were subjected to acidolysis conditions without EG as a stabilizer (Figure , path 1); in the second case, EG was used as a stabilizer (Figure , path 2). A typical reaction was carried out as follows: substrate (1-G or 1-S, 1 equiv, 0.024 mmol, final concentration in DMC: 0.016 mmol/mL), EG (0 or 4 equiv, 0 or 0.096 mmol, path 1 and 2, respectively), internal standard (1,2,4,5-tetramethylbenzene, 0.04 mmol, 1.7 equiv), DMC as solvent (1.5 mL), and H2SO4 (5 or 2 mol % to the substrate, from a stock solution in DMC of concentration 0.1 mmol·mL–1) were added to a MW vial. The vials were sealed, purged with nitrogen, and the solutions stirred at 650 rpm and heated to 140 °C in a heating block equipped with an extra vial with a thermometer as control for the needed time. At the end of the reaction, the vial was cooled down in an ice-bath, and a sample of 0.35 mL filtered through Celite and NaHCO3, and either injected directly into the GC-FID or diluted in 1.3 mL CH3CN for HPLC analysis (see SI).

Computational Methodologies

The computational procedure consisted of (1) running atomistic molecular dynamics simulations based on a reactive force field (ReaxFF)[67] to sample possible configurations of the model compounds in a DMC solution containing catalyst and stabilizer (Figure S1 and SI); (2) extracting the most promising structures where at least one molecule of both the catalyst and stabilizer were within 3.5 Å of the OH moiety of the solute; (3) using these configurations, appropriately reduced, as starting points of quantum chemistry (QC) density functional theory (DFT) calculations. This last step had the purpose of identifying possible reaction pathways through the nudged elastic band (NEB) methodology, implemented in the Quantum Espresso package (QE),[68] and give an estimation of the activation energy barriers. All the details are reported in SI.

Possible reactant and product geometries extracted from the simulations were optimized at the DFT level and used as starting and final structures of the NEB reaction paths. The minimum energy paths (MEP) between the two minima and the saddle geometries were identified using the NEB method with convergence criteria of 0.1 eV Å–1, considering 25 intermediate images between the local minima in each case. All DFT calculations used plane waves UltraSoft (US) pseudopotentials and the Perdew–Burke–Ernzerhof (PBE) functional. The energy cutoffs for the wave function and charge density and potential were 40 and 400 Ry, respectively. Gaussian smearing of about 0.002 Ry was used in all the calculations that were spin-unpolarized. Brillouin zone sampling was restricted to the gamma point. The Grimme-D2 method, which accounts for dispersion forces, was employed. Other details are reported in SI.