Modelling the fragmentation kinetics of the heterogeneous lignin macromolecule during kraft pulping with stochastic graphs

This article presents a novel concept for modelling the kinetics and related phenomena of the kraft pulping process on a macromolecule level with the initial objective to explicitly model and relate the breakage of phenolic and non-phenolic β – O – 4 bonds to the observed three-stage delignification profile. The modelling frameworks consist of the building and the fragmentation of the lignin macromolecules. The macromolecules are modelled as stochastic graphs where monolignol object nodes are reassembled in a Monte Carlo approach into internal structures, which aggregate to a lignin macromolecule interconnected through chemical bonds represented by the edges of the graphs. The fractionation follows the splitting of β – O – 4 ether bonds with different configurations resulting from functional groups attached to the monolignols, namely the phenolic and non-phenolic β – O – 4 bonds with their respective stereochemistry. It is tested against a previously published model based on an extension of the established Purdue kinetic model and experimental data. The results align with the observed delignification trajectory during kraft pulping, and the hypothesis that β – O – 4 bonds splitting is mainly responsible for the delignification. However, some discrepancies between the current model, the previous model and experimental data are presented. These differences are discussed in the context of recent experimental findings indicating that β – O – 4 bonds splitting might not entirely be responsible for the delignification due to mass transfer


Introduction and motivation
With the acceleration towards an economy based on green and renewable feedstocks, attempts have been made to maximise the utilisation of lignocellulosic raw materials.In particular, the valorisation of lignin poses one significant challenge in today's research.With this transition, modelling approaches can help guide biorefinery processes to optimise their utility.However, lignocellulosic feedstocks are inherently complex on several levels, and their understanding is crucial for model development.
On the macromolecule level, lignin is often described as a highly heterogeneous macromolecule with its three main phenylpropane monomers (monolignol), paracoumaryl alcohol (H-unit), coniferyl alcohol (G-unit) and sinapyl alcohol (S-unit), in different proportions depending on its feedstock origin.These monolignols are linked via biosynthesis through various chemical bonds such as ether bonds: β -O -4 (most abundant bond) etc. or condense structures such as 5-5 bonds, etc., to internal structures (see Fig. 1), which in turn are aggregated to a larger macromolecule (see Fig. 2).The monolignols within the macromolecule can come in different chemical configurations regarding their functional groups, e.g.phenolic or non-phenolic monolignols.Moreover, lignin is considered to be linked to carbohydrates through lignin-carbohydrate complexes (LCCs) presented within the lignocellulosic material.[1].
The current structural analysis is mainly based on wet chemistry and NMR analytics of Björkman lignin [3], also called milled wood lignin (MWL).The MWL is produced by mechanical milling and extraction via an aqueous dioxane solvent.The obtained isolated lignin is considered the closest approximation of lignin compared to its native state.However, the ball milling does affect the lignin structure, e.g., a change in β-O-4 bond content [4] and molecular size [5].
New fractionation protocols have been developed to dive deeper into the lignin heterogeneity, addressing questions regarding possible changes in native lignin during ball milling.Essential remarks regarding the origin of MWL, whether the lignin originated from the middle lamella [5,6] or the secondary wall [7], the association between lignin and carbohydrates, and whether lignin is branched [2,8] or a linear [9][10][11] macromolecule have been made.[8] Given the inherent heterogeneous nature of lignin, another complication occurs during its fractionation.Chemical fractionation processes like kraft pulping alter the chemical structure of native lignin severely [12].In terms of lignin valorisation, these fractionation processes represent the first stage of the lignin supply chain, followed by possible depolymerisation and upgrading of the fractionated lignin [13].They play a crucial role in the fate of the lignin valorisations and need to be considered carefully.
Recent experimental research has shown that lignin from different sources/fractionation processes experiences molecular weight and functionality differences, which correlate with differences in, e.g.mechanical, electrical, or thermal properties.Hence, attention should be paid to the structural and chemical properties of the lignin during the engineering of lignin-based materials [14][15][16][17] and theoretical work on modelling the lignocellulosic feedstock and its fractionation must adjust to the experimental resolution [18].
A primary source of technical lignin comes from the pulp and paper industry in the form of kraft lignin.Advancing and optimising pulp mills with an additional objective to produce well-defined lignin fractions starting from the fractionation process over subsequential downstream processes would improve their productivity and contribute to the ambition of a sustainable future.[19].
Reviewing the theoretical approaches to predict the fractionation of lignin during the kraft pulping process leaves one with two prominent kinetic models, which treat lignin as a homogeneous chemical component.The first model considers the fractionation process in a three-stage delignification model, and the second includes high and low reactive lignin portions [20][21][22][23][24].Both models disregard the lignin macromolecule's heterogeneous nature, lumping its complexity into a small number of parameters while implicitly assuming differences in lignin reactivity due to its chemical structure.Differences in splitting kinetics of the lignin interunit linkages and solubility/mass transport phenomena of different molecular sizes are entangled and cannot be modelled explicitly with the widely used approaches.Newer modelling approaches have tried to include the heterogeneous nature of the lignin to some extent, where the delignification is considered with a continuous reactivity distribution type model [25] or a more detailed description for the fractionation of lignin in a microkinetic approach [26].These models also concentrate on the delignification and the removal of lignin from the feedstock rather than enabling the prediction of lignin fragments on a macromolecule level with an entire decoupling of kinetic, solubility and mass transfer effects.
Recent multiscale modelling approaches have tried to include some heterogeneous effects during kraft pulping on a cell wall/ fibre-to-fibre resolution using the Purdue kinetic based on a Monte Carlo kinetic approach, considering differences in concentrations of the respective morphological region [27][28][29].However, the heterogeneous nature of lignin on a macromolecule level is still omitted.
This article provides a model concept for the case of softwood kraft pulping delignification to build the bridge from the current lump kinetic approaches to a more systematic mechanistic model, which enables the prediction of the heterogeneous nature of the lignin fractions on a monolignols resolution, shifting away from the lump kinetic approach, as mentioned in earlier publications and researchers [30][31][32].The initial objective is to model the lignin fragmentation by explicitly relating the macromolecule's fragmentation to the splitting of phenolic and nonphenolic β -O -4 bonds, considered the dominant scissoring reaction  during kraft pulping.Moreover, it aims to disentangle the kinetics of bond breakage from other phenomena, such as mass transfer and solubility, which have lately been hypothesised to influence the delignification [33][34][35][36][37][38][39].

Building the lignin macromolecule
The first step of developing a kinetic model that describes the kraft pulping fractionation process on a monolignols resolution is to define the lignin macromolecule.Earlier researchers attempted to model lignin biosynthesis in silico with different approaches, from combining Monte Carlo and first principal methods to force field computation [40][41][42][43].This article will focus on a more practical approach to generate the lignin macromolecule based on a model for spruce MWL presented by Balakshin et al. [2].The model considers lignin as interconnected monomer units.It follows the connection of these units via chemical bonds, which are represented as edges in a graph theory approach.Lignin growth is modelled by appending monomers using a Monte Carlo approach where the possible connections between monomers are predefined by their tail and head portions.Afterwards, these object nodes are stored in a graph structure, where the edges represent chemical bonds, and the nodes are monomer objects with, e.g., different functional groups at a given location on the monomer.Since we will consider softwood as the lignin source, only coniferyl alcohols (G-unit) are currently considered in the model, neglecting the small percentage of paracoumaryl alcohol (H-unit) presented in softwood to simplify the model.Fig. 3 illustrates one G-Unit as an object used in the model.
To the commonly recognised lignin structures (Fig. 1), a subset of non-resolved alkyl-O-alkyl ether structures, considered to give rise to branching, are included to be consistent with the Balakshin et al.MWL model [2,44].The conceptual algorithm for the lignin macromolecule generation is given in Fig. 4.
Initially, the model decided from a Log-normal distribution the size of the macromolecule.Afterwards, the first lignin monomer is drawn from the monomer library according to its probability and added as a node to a graph representation.Monomers with open connections within the graph are tracked.A random node from the monomers with an open connection is chosen, and a fitting counterpart is added according to its possible connection property.Head-to-head, tail-to-tail or head-to-tail connections determine the possible connections.Fig. 5

offers two examples of possible connections between monomer objects.
Monomers are added to the macromolecule until the predefined macromolecule size is reached.The Supplementary material provides more details with an example of the lignin macromolecule generation and all monomers used in the model with a description of the nomenclature (Supporting materials A and C).

Fragmentation of the lignin macromolecule
The second step in the model is the fractionations of the created lignin macromolecule due to the splitting of chemical bonds.Fig. 6 presents the conceptual framework for the model on three different levels.The lowest resolution of the model is represented by a softwood fibre with its three distinct phases: fibre wall (solids and fibre wall liquor), fibre lumen and free liquor.The model's frame will be the delignification within an element of the fibre wall.Zooming into this fibre wall element (white box in the schema), the lignin macromolecule reassembles monolignol building blocks connected through chemical bonds (edges in the graph representation).On the third level, these blocks are represented by their chemical bonding and functional groups (object nodes in the graph representation).
The model aims to initiate the disentanglement of kinetics and solubility effects, shifting away from lump kinetic approaches.During the fractionation, individual chemical bonds and their breakage are followed.The model's initial objective is to incorporate studied splitting reactions of the β-O-4 ether bond, considered the dominant scissoring reaction.To include some nuances of the β -O -4 splitting, the following four configurations of the ether bond are included in the model, phenolic and non-phenolic β -O -4 bonds with their respective stereochemical configurations erythro and threo, presented in equal amounts.The kinetic parameters (Table 1) for the splitting are taken from the experimental investigation by Ahvazi and Argyopoulps [45].The phenolic β -O -4 bonds follow the parameters for the splitting in the initial delignification phase, and the non-phenolic β -O -4 bonds align with the bulk delignification parameters of the study.The solubility of a lignin macromolecule results from its breaking up into smaller lignin fragments.Fig. 7 illustrates schematically how the delignification in the model progresses.
Time step t 1 represents the initial solid lignin macromolecule within the fibre wall.During a time interval, Δt an arbitrary number of bonds can break, where the breaking events are independent of one another.All fragments below a predefined threshold macromolecule size are dissolved.The model progresses in Δt time steps until the end time is reached.
The splitting of each β-O-4 bond configuration is assumed to follow a pseudo-first-order reaction with respect to the chemical bond.
The probability of breaking a bond p β−O−4i is expressed by the conversion of the bond type: Here Δt is the time step and k i is the reaction constant calculated according to the Arrhenius equation: The yield of one lignin macromolecule is defined as: It is the quotient of the non-dissolved nodes and the initial presented number of nodes in the macromolecule.Further assumptions used in the model are: N. Bijok and V. Alopaeus • Only molecular size dependency for the solubility is considered.However, other effects, such as pH or temperature, could be incorporated into the model with rigorous conducted experimental data [46,47].Experimental research gives a range for the number average molecular weight of different kraft lignin fractions between M n = 900---2000 g/mol with an M n of 1400 g/mol for the unfractionated kraft lignin macromolecules [12].In another article, the M n is given as 1587 g/mol [46].The current model uses a 1440 g/mol M n , equivalent to an eight-node lignin macromolecule.• In the first approximation, the model assumes that the β -O -4 bond cleavage is the leading cause of delignification, omitting the breaking of the other chemical bonds.• After the lignin is soluble, the mass transfer of the macromolecule from the fibre wall liquor to the free liquor phase occurs immediately, assuming no mass transfer limitations.• The model uses the MWL model suggested by Balakshin et al. [2] as an approximation to represent lignin's native state.
• Only lignin is considered in the model, and possible LCCs are not presented in the macromolecule.• The kinetic parameters from Ahvazi and Argyopoulps [32] are investigated at an effective alkali concentration of 15 % on wood and a sulfidity of 13 % in a liquid-to-wood ratio of 4/1.Hence, the pseudo-first-order reaction rate constants are valid for these kraft pulping conditions.

Results and discussion
According to the algorithm presented in Fig. 4, lignin macromolecules have been generated, and the lignin structure properties are presented in Table 2 as structures per 100 C 9 -units.An optimisation has been conducted with the particle swarm optimisation algorithm from the MATLAB optimisation toolbox to match the experimental data from Balakshin et al. [2].However, no attempts have been made to match the experimental data perfectly.As mentioned by the authors, these lignin  N. Bijok and V. Alopaeus structures represent a snapshot of lignin molecules within a broader distribution of possible lignin macromolecules.Hence, the given data/ model does not necessarily represent the average lignin molecule but produces reasonable macromolecules that match overall experimental observations.The Supporting materials B describes the optimisation procedure.Fig. 8 illustrates an exemplary eight-node fragment of the lignin macromolecule from the graph representation used in the current model.
The average lignin yield profile of simulating several lignin macromolecule realisations (solid line) and their relation to the respective β-O-4 bonds splitting (see Fig. 10) is compared to the average delignification profile of our previous published model [31,32] with the Ahvazi and Argyopoulps [32] kraft pulping condition (dot-dash line) for small wood chips (5mmx2mmx2mm) and experimental data from Aurell and Hartler with similar effective alkali charge but higher sulfidity (effective alkali 15.75 per cent on wood, sulfidity of 25 per cent and a liquid-to-wood ratio of 4/1) (dots), showing the typical three-stage delignification profile during kraft pulping [48] (see Fig. 9).
The lignin yield over time for the current model follows the observed S-shape curve during kraft pulping.In the initial phase, the delignification ramps to around 100 min, corresponding to the splitting of mainly phenolic β -O -4 bonds (see Fig. 10).The bulk delignification starts at around 100 min (delignification of 20 %) at around 150 • C. A sharp decline in lignin can be seen in the bulk phase, which starts to slow down when the delignification reaches about 50 per cent.At about 200 min, the delignification enters the residual phase (delignification of 90 %).At the end of the simulated time, most of the lignin is delignified (around 95 %), and a small percentage of the non-phenolic β-O-4 bonds (around 4 %) haven't reacted.
Comparing the current model with our previous modelling approach Fig. 6.Framework of the model including three levels.First level: Softwood fibre divided into three phases, fibre wall, fibre lumen and free liquor.Second level: The lignin macromolecules build up of monolignol building blocks.Third level: Individual monolignols with their chemical bonding and functionality.

Table 1
Kinetic parameters for the β -O -4 bond splitting according to Ahvazi and Argyopoulps [45].N. Bijok and V. Alopaeus [31,32], where we used the Purdue kinetic model as an initial starting point and assigned an empirical yield-dependent reactivity decline factor, which assumed a change in lignin reactivity on a chemical structure level, yields a similar delignification trajectory as in the current model.Both models align almost identically in the initial and first part of the bulk delignification phase, with a noticeable deviation between the models at around 50 per cent delignification and reaching the intended cooking temperature.From the second part of the bulk delignification, the previous model approach converges to a higher residual lignin percentage than the current model.When comparing both models to the experimental data [48] with similar effective alkali concentrations but higher sulfidity, the deviation prevails in the initial part of the delignification.This is in line with the current understanding of the reaction mechanism of phenolic β -O -4 bonds, which splits in the presence of HS − ions [49].On the other hand, the experimental data behaves similarly to our previous model approach in the residual phase, yielding a higher residual lignin content at the end of the cooking process compared to the current model.
Overall, the current model predicts well the generally observed characteristics of the kraft pulping delignification in a mechanistic approach by explicitly relating the delignification trajectory with the different β -O -4 bonds configuration cleavage instead of using a yielddependent reactivity decline factor as previously suggested [31,32], implicitly assuming the different β -O -4 bonds splitting with a three stages kinetic model [24] or lumping the lignin into low and high reactive parts [20][21][22][23].However, there are some discrepancies in the residual lignin content in the later part of the delignification compared to the experimental data and the previous models.
Lately, the hypothesis that the β-O-4 bond splitting is entirely responsible for the delignification has been challenged.Moreover, it is assumed that some β-O-4 bonds are not reactive towards the kraft pulping conditions and/or mass transport/solubility phenomena affect the delignification process [33].Experimental research shows that different lignin macromolecule fragments have mass transport limitations through the fibre cell wall.This results in mass transport limitations from the fibre cell wall to the free liquor presented in the digester [34][35][36].A recent study showed that lignin macromolecule fragments could be leached out from unbleached pulp before oxygen delignification, indicating that soluble lignin macromolecules are presented in the cell fibre, but mass transport phenomena prohibit the delignification of these fragments during the pulping process [39].Additionally, LCCs linked to the lignin macromolecule are assumed to be presented in wood chips, contributing to delignification resistance in the residual phase of the process [50,51].
The current model is limited to the hypothesis that the β -O -4 bond splitting is entirely responsible for the delignification.However, the framework allows adding a sub-model for mass transfer effects on a fibre cell level where lignin macromolecule fragments of different molecular sizes can be assigned with varying mass transport properties [34][35][36].Related to that, the solubility properties of lignin, e.g.due to different N. Bijok and V. Alopaeus ionic strengths impacting the molecular weight distribution of the lignin fragments, can also be considered with additional experimental investigation [38].Extending the graph representation of the lignin with LCCs monomer objects connected to the monolignol and their influences on the residual lignin can be explored with the model structure [8].Finally, adding reaction mechanisms on a functional group/chemical bond level, e.g. the breaking of other interunit bonds, the formation of enol ether structures leading to alkali-stable structure [52] and possible condensation reactions [53,54] acting on the monomer objects will help to predict the resulting lignin fragments in more detail, especially in the later bulk and residual phase of the process.
With the given resolution of the model, an interesting comparison between experimental research and model results can be drawn regarding the topochemical behaviour of lignin.Fig. 11 compares the β-O-4 bond fractions over the delignification degree of the different β -O -4 bond configurations with topochemical experimental results from Saka et al. [55] for the lignin in the cell corner (CC), compound middle lamella (CML) and secondary wall (S).As shown in Fig. 11, up to 20 % delignification, the phenolic β-O-4 bonds are primarily split, corresponding in the topochemical experiments to a slight delignification in the secondary wall.When the delignification reaches around 20 % (end of initial delignification phase), the non-phenolic β -O -4 bond starts to split.This corresponds to the delignification in the cell corner and compound middle lamella.The authors speculated that the secondary wall is higher in phenolic β-O-4 units, whereas the CC and CML contain more non-phenolic β -O -4 units.This has been observed   Given the discussed limitations, further exploration of the topochemical effects, reaction mechanism and mass transfer/solubility effects will help to guide a more target-oriented delignification of specific regions and defined lignin fragments from the lignocellulosic feedstock.The ultimate goal is to predict lignin fragments with their difference in the functional group, yielding a targeting use of the lignin fragments in the next steps of the valorisation supply chain [14][15][16][17].

Conclusions
This article presented a new concept for modelling lignin fragmentation during kraft pulping using stochastic graphs, which enhances the possibility of drawing conclusions from the kinetic data based on welldefined compounds to real delignification during kraft pulping and other relevant biorefinery processes.It's shifting away from lump kinetic towards an approach that allows a deeper understanding of the individual reactions in delignification and enables us to model the delignification processes more systematically and mechanistically by considering lignin's heterogeneous nature.The initial objective of the model was to enable the decoupling of reaction kinetics and solubility phenomenon, allowing the prediction of the delignification on a monolignol resolution.The model's current state has been limited to splitting the β-O-4 bond of phenolic and non-phenolic lignin structures.Moreover, simplification regarding solubility has been made.Nonetheless, the kraft pulping process simulation with the described assumptions agreed qualitatively with the commonly observed experimental delignification trajectory.
However, additional reaction mechanism that alters the lignin fragments and the mass transfer of these fragments from the cell wall to the free liquor, which are hypotheses to also play a role in the overall delignification, have been omitted at this stage and must be further explored to include these effects into the model.
In summary, the article presents an open-ended kinetic model approach providing a theoretical tool to decipher the lignin puzzle and its fragmentation during kraft pulping.The limitation of the model has been pointed out and discussed based on experimental research, branching into further question that needs to be addressed and investigated to generate a more comprehensive fractionation model.Filling in the discussed limitation and keeping the model updated with new experimental findings will facilitate the ambition to predict and produce well-defined lignin fragments from kraft pulp mills.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 5 .
Fig. 5. Examples of monomer connections.The head-to-tail connection between monomer 1 and 2 via β -O -4 internal structure and the head-to-tail connection between 2 and 3 via Phenylcoumaran (left side).Head-to-head connection of monomer 1 and 2 via 5-5 bond and head-to-tail connections of 3 with 1 and 2 via β-O-4 and α-O-4 bonds to form dibenzodioxocin (DBDO) internal structure (right side).

Fig. 7 .
Fig. 7. Schematic illustration of lignin fragments' chemical bonding splitting and solubility for four different time steps (t 1 -t 4 ).Encircled graphs represent dissolved lignin fragments in the next time step.

Fig. 8 .
Fig. 8. Illustration of the graph representation for a lignin molecule fragment containing eight lignin monomer units (node six still has an open connection for a 4-O-β unit, and node eight can add a β-O-4 and 5-5 unit.Additionally, nodes 6 and 8 lack a DBDO/α unit to form a cyclic dibenzodioxocin (DBDO) internal structure).

Fig. 10 .
Fig. 10.Bond splitting profile over time for the different β-O-4 bonds included in the model.

Fig. 11 .
Fig. 11.Comparison between the bond splitting simulated in the current model and topochemical delignification behaviour of lignin in cell corner (CC), compound middle lamella (CML) and secondary wall (S) [55].

N
.Bijok and V. Alopaeus

Table 2
[2]uctural description of the lignin macromolecule according to Balakshin et al.[2]and the current model.