Experimental validation of a reaction network model for autoxidation of linoleate esters

A deterministic mathematical model is used to study complex chemistry and determine kinetic parameters of the autoxidation polymerization (drying) of fatty acid esters such as ethyl and methyl linoleate (EL and ML). The model was developed previously and is based on an algorithm that generates a reaction network and the associated kinetic equations in an automated way (automatic reaction network generation – ARNG). The ARNG model computes the concentrations of the monomeric species, including the number and type of crosslinks (peroxyl, ether, alkyl) they possess. In this study, kinetic parameters of the main reaction families were estimated by adjusting them to obtain a good fit with new experimental data of various nature. Among these were: Fourier transform infrared (FTIR) spectra, electron spray ionization mass spectrometry (ESI-MS) data of EL with TiO 2 cured at 70 ◦ C and size exclusion chromatography (SEC) data of ML cured at 80 ◦ C. Sensitivity analysis revealed that addition, bis-allylic hydrogen abstraction by peroxyl radical and peroxyl radical recombination are the most influential reactions in the considered set of functional groups. These sensitivity results indicated which reaction families were considered in the parameter estimation procedure. Rate coefficients of key reactions like hydrogen abstraction, O 2 to radical addition, OO-and O-radical addition to conjugated double bonds and β -scission could be well be determined from EL-conversion and conjugated double bond profiles from FTIR. Some tens of dominant peaks in the ESI-MS spectrum were also found in the ARNG model, while vice versa most dominant model predicted peaks were confirmed by ESI-MS. ARNG allowed detailed chemical pathway analysis in terms of intermediate products, which was helpful in elucidating the pathway for the most important acidic monomer, hexanoic acid. Finally, the most dominant crosslink type formed turned out to be peroxyl, less ether, while alkyl is practically absent. The ARNG model of linoleate drying has thus considerably gained power to predict properties relevant to the chemical structure of EL polymer.


Introduction
Natural oils are complex chemical mixtures.For many years, the drying oils that comprise the binder portion of paints have been studied under ambient conditions using a variety of approaches, both experimental and computational [1][2][3][4][5][6][7][8][9][10].Several factors drive the motivation of these studies, like for instance in the area of oil painting conservation, where they contribute to understanding chemical and physical mechanisms of paint degradation.Other emerging applications of natural oils are found in food technology [11][12][13][14]and in development of renewable feedstock for polymers, replacing fossil fuel resources [14][15][16][17].In this context, linoleic acid peroxide obtained by autoxidation of linoleic acid is used as macro-initiator for common polymers [15].Schaich [5,6] and Hazer [14] provide reviews from a food technology perspective, the latter being concerned with the use of linseed oil autoxidation as a marker to test the antioxidant effect of packaging polymers like polypropylene grafted with substances like tannic acids.Orlova's [18] review covers the oil paint area, both from experimental as well as modeling perspective.
Drying oils, such as linseed oil, cure via autoxidation, a free-radical polymerization process through which oxygen adds to the oils after an initiation event [5,6].Subsequent oxidative reactions occur that diversify the distribution of products alongside crosslinking reactions between the alkyl tails of triglycerides that lead to dense, polymer networks [19][20][21][22].Specific observable phenomena are a consequence of the autoxidation chemistry in these oil systems, such as smell during drying due to volatile organic compound formation [8] or visual changes during ageing due to metal soap formation.In recent years the latter phenomenon has become of interest, ever since it has been shown that there is a strong link between the formation of these metal-fatty acid complexes and the concentration of carboxylic acid groups.In ageing oil paintings acid concentration further increases due to extended oxidation and also by hydrolysis of the ester bonds [23].Thus, having a detailed sense of how key functional groups react in a given system can improve understanding of long-term behavior of oil paints, leading to more optimal conservation strategies for art objects.
Unlike the regular constituent monomers in many industrial copolymers like ethylene-propylene-third monomer (EPT) or acrylonitrilebutadiene-styrene (ABS), the polymer network of dried linseed oil (LO) consists of connected monomeric units of a wide variety, produced by the autoxidation process, indeed much wider than industrial copolymers.Although linseed oil homopolymerization starts off with just a single LO monomer, each molecule has multiple reactive sites of varying propensity on each of its three fatty tails.During polymerization the autoxidation leads to an explosion of different monomer units formed over time, which crosslink in innumerable ways to form the dried polymer network.Simplified systems based on pure fatty acid methyl/ ethyl esters have also been studied that capture representative reactivity without the variation in degree of unsaturation as found in drying oils [24][25][26].Fatty acid esters are reasonable model molecules for synthetic alkyd coatings, as they form mainly oligomers and extend the investigation of the curing of coatings beyond traditional, bio-based oil origins [27][28][29][30].The complex mechanisms of the underlying autoxidation process of drying oil have been studied by numerous authors for decades, most comprehensively by Schaich [5,6], but these remain qualitative.Recently, quantitative kinetic models became available, mostly from our group [9,18,31,32], but these were not or only partly validated with experimental data.In the present paper we report a considerably better validation of a recently published model based on the ARNG methodology [32].
The main purpose of this work is to develop a comprehensive and validated model that describes reactions involved in drying of oil paint.It is based on recent advancements in computational modeling using a monomer approach with automated reaction network generation (ARNG) applied to the curing of complex drying oils by Orlova and colleagues [32].At its foundation, molecular graphs representing monomers are modified based on reaction rules until no more unique monomer units, both connected and non-connected, are found.Due to the complexity of the autoxidation process, even under simplifying assumptions imposed on the reaction rules, thousands of monomer units in a reaction network comprising millions of reactions are found.The reaction rules represent the autoxidation chemistry of drying oils without catalyst and cover most of the mechanisms known in literature, most comprehensively described by Schaich [5,6].In the ARNG model the concentrations for the mentioned monomer species or their functional groups are solved with a continuum kinetic model using thousands of kinetic parameters.After the monomer concentrations are solved, polymeric species can be assembled using random graph (RG) theory to connect all the "half-bonds" on monomers until it is a fully converged system.This is an effective method to uncover certain properties of highly crosslinked polymers.However, the original ARNG model was not tuned to experimental data.In fact, it is virtually impossible to validate all the thousands of parameters of the ARNG model.In this work we have employed both previously published data on the drying of EL as well as new data generated for the purpose of this publication.Data on EL autoxidation is rather rare, since uncatalyzed drying of linoleate esters is very slow at room temperature, with moderate acceleration at elevated temperatures.We will refer to studies of linoleic acid and EL autoxidation by Hazer [14], Holman and Elmer [33], Oyman et al. [34][35][36], Muizebelt et al. [37,38] and a size-exclusion chromatography (SEC) analysis of ML by Mourão [39].Our own new experimental data on concentrations of certain key components, like the conversion of the starting monomer EL, were generated by Fouriertransform infrared spectroscopy.Together with the older data we have identified kinetic parameters that dominate the key component's concentration trends.
The estimation of kinetic parameters in the ARNG model faces several challenges, which are partly inherent to all larger parameter estimation problems.On the one hand, the complex reaction network has many parameters, possible too many in view of the limited set of experimental data available.On the other hand, the found estimates might contradict theoretically expected values, which then would ask for additional explanation and possibly extensions of the already large model.The parameter estimation in this work has been executed according to a strategy, where first various sets of included chemical reactions were tested, by varying the number of reaction rules applied and observing the impact of that on the experimentally verifiable model outcomes.Subsequently, the rate parameters for the optimal set of reactions were determined by varying their values and comparing model predictions with experimental data.A summary of experimental concentration profiles is given, predominantly from attenuated total reflection Fourier transform infrared (ATR-FTIR) spectroscopy.Furthermore, we perform an analysis of the sensitivity of functional group concentrations to rate constants of key reaction steps and list final optimized kinetic parameters for lumped H-abstraction rates, hydroperoxide decomposition, and radical termination.Finally, an analysis of electrospray ionization-mass spectrometry (ESI-MS) mass spectra peaks is performed in both positive and negative modes.We find an overlap between model-generated monomers and dimers and experimental data, and the highlight regions of measured data where the model still lacks precision in its chemical space.

The ARNG model
The kinetic model follows the framework by Orlova et al. [32].Mathematically, molecules are represented as graphs with atoms as nodes and covalent bonds as edges.In this representation, functional groups (e.g.alcohols, hydroperoxides) are subgraphs or certain substructures of larger molecular graphs.These functional groups, or reactive patterns, as they are called by Orlova et al. [32], for autoxidation chemistry are shown in Supplemental Fig. A.1.The subsequent model is also used in this work.Molecular graph connectivity is stored in an adjacency matrix that weights bonds by their bond order.Radicals, as unpaired electrons, are denoted on the diagonal of the square adjacency matrix.
For polymerization, it becomes increasingly infeasible to track the connectivity of each atom in the growing polymer.Thus, Orlova et al. [32] subdivides the polymer into monomeric species.Rather than crosslinking into larger polymers, the species only maintain the position of a crosslink and its type, in addition to other local functional groups.For the case of drying oil autoxidation, 3 crosslink types are considered: alkyl, ether and peroxyl.One advancement made to this model since publishing is the directed crosslinks for the monomers to distinguish crosslinks made by radical-radical termination reactions (non-directed crosslinks) and radical-to-conjugated double bond addition reactions (directed crosslinks), depicted in Fig. 2.This distinction may influence the global polymer properties that are inferred from the crosslink distribution using the RG-model [40][41][42] indeed, the idea of crosslink directionality is adopted from directionality of edges in graph theory.Although modeling the polymer network formation is mostly beyond the scope of this work, we did use the RG-model to investigate the sensitivity of crosslinks and polymer network gel formation to certain rate coefficients (Paragraph 3.5).
The monomeric species are generated on-the-fly using the ARNG model, a flowchart is shown in Fig. 1.Reactive patterns (see Fig. A1, SI) are identified in the molecular graph (see Fig. 1a) and are transformed by adding or removing edges, equivalent to the making and the breaking of chemical bonds.This happens according to reaction rules or set transformations being provided to the model a priori.In an automated fashion, the reaction rules are applied to each reactive pattern of a reactant molecular graph, which is transformed into a product molecular graph.Then the algorithm checks if such product has already been formed.This process iteratively fills in the potential chemical space.Thus, the ARNG model deals with the reactive parts of the EL and ML molecules.All these transformations are stored in a reaction network that captures the reactants, reactions, and products.This network is then transformed into a system of ordinary differential equations, where each equation describes the change in concentration for each monomeric species.The rates of change in concentration in these equations are parametrized by rate coefficients that are either experimentally found or derived via the Arrhenius equation.The concentrations of each monomeric species and each functional group of interest are calculated and lumped for comparison with experimental results.To decrease computational complexity of the algorithm, the graph representation in the ARNG model omits the ester head and the carbon tail (see Fig. A2 of SI) that do not take part in the autoxidation process.When analyzing molar masses of the generated monomeric species, head and tail masses are added back in to end up with the real molecular masses.
Concerning rate coefficients, a simplification with respect to the approach of Orlova et al. [32] has been applied.Originally, all coefficients were calculated using group contribution theory, which gave rise to huge variations within one reaction family, often even by more than ten decades.We consider these variations as unrealistic and, abandoning this approach, instead have employed average rate coefficients for each reaction family.The average values still find their frame of reference in structure-reactivity relationships, as we will discuss in 'Results'.Thus, the coefficient values displayed in Table 1 are the values used in the ARNG model and are no longer modified by group contribution parameters.Also, as a consequence, we were able to strongly simplify the peroxide crosslink decomposition step, since the decomposition of a crosslink of any monomeric unit can be described as a first-order reaction just dependent on the concentration of that unit and a constant rate coefficient.This rate coefficient is no longer dependent on its structure nor on the structure of the connected unit, which leads to considerable reduction of computation time.

Mechanisms
The autoxidation reactions for a free-radical polymerization included in this model are summarized in Table 1 of the Results section.Beyond these specific reactions, there are reactions of secondary products where reactive patterns are identified to account for greater complexity.The model is extended as compared to Orlova's model [32] by including reactions with the highly reactive hydroxyl radical: radical addition, H-abstraction from bis-allylic position and HO• termination with itself [43].As in Orlova et al. [32] it is assumed that crosslinking takes place by radical-radical termination and addition reactions of peroxy-and alkoxy radicals to conjugated double bonds.This implies  The two acids have the same reaction product part (species nr1009), but a different tail.Type of reaction is shown along with total amount reacted from reactant(s) to product(s) in 25 h.Fig. A16 of SI shows all the species' concentration profiles over time.EL (species 3), original concentration 3.71 mol/L, fully reacts away by Habstraction(HA) to radical species 5.After a number of reaction steps aldehyde species 531 is produced by β-scission type 2. Just over 1 % of species 531, 1.48⋅10 − 5 mol/L, finally produces the acids through Baeyer-Villiger reactions with all peroxyl acids.The most important acid within this group of peroxyl acids is species 1963, responsible for 13 % of the production of hexanoic and 9-ethoxy-9-oxononanoic acid.

Table 1
Autoxidation reaction families and their respective rate coefficients at 70 • C in three sets.Values with * optimized in this study.

Nr
Reaction family Functional group representation Secondary initiation by initiator radical IOO  that three types of crosslinks are possible: alkyl, ether, and peroxyl.Note that termination and radical addition reactions not involving HO• Radicals contribute to polymerization and network formation, while HO• radicals do not.

Experimental techniques 2.3.1. ATR-FTIR
An experimental study was conducted with a mixture of ethyl linoleate (Sigma-Aldrich) and coated rutile titanium dioxide pigment (Tronox CR-826).This pigment is considered to be chemically inactive during drying and degradation [44] ATR-FTIR spectra were collected on a Frontier spectrometer (PerkinElmer) equipped with a heated diamond GladiATR module (PIKE Technologies).A drop of EL was placed on the ATR-diamond and subsequently heated at 70 • C during spectrum collection under normal lab atmosphere (50 % relative humidity).Spectra were collected every 12 s at 4 cm − 1 resolution.Concentration profiles of several functional groups were obtained by extracting the FTIR intensities of their corresponding absorption bands after baseline correction: carboxylic acid (1715 cm − 1 , 1412 cm − 1 ), conjugated double bonds (987 cm − 1 ), single C --C bonds (975 cm − 1 ) and a group around 3100-3600 cm − 1 , attributed to alcohol, carboxylic acid and hydroperoxyl (Fig. A3 in SI).

ESI-MS
The same system of ethyl linoleate and coated rutile TiO 2 was used for conducting mass spectrometry experiments.The paint was kept at room temperature in lab conditions and then at 70 • C for 146 h without light, and samples were taken at 0, 24, and 146 h.The experiment was performed from start to finish three times and the same trends in results were obtained, taking consideration that exact sample times could not be repeated as it is not an online measurement.
The extractable fraction, mainly of monomers and dimers, was studied by direct-injection ESI-MS measurements taken in both positive and negative modes with a mass resolution of 0.15 m/z.Experimental details and methods for peak analysis follow those published in the Orlova et al. [32] study on the autoxidation of triolein [45].

Results
The main objective of this work is to validate the previously developed model of oil paint with experimental data.This is a great challenge as the ARNG model contains a great number of kinetic parameters to be estimated for the reactions between thousands of species.For validation purposes, we have employed both previously published and newly generated experimental data on the drying of EL.We refer to studies of EL and linoleic acid autoxidation by Hazer et al. [15][16][17], Holman and Elmer [33], Oyman et al. [34][35][36], Muizebelt et al. [37,38] and a sizeexclusion chromatography (SEC) analysis of ML by Mourão [39].Our own new experimental data on concentrations of certain key components, like the conversion of the starting monomer EL, were generated by Fourier-transform infrared spectroscopy.Experiments by Hazer et al. [15][16][17] revealed autoxidation of linoleic acid at room temperature, while both EL and linoleic acid were shown to undergo autoxidation at 37 • C by Holman and Elmer [ 33].Latter authors observed that autoxidation rate of acids is significantly higher.The EL experiments by the Oyman and Muizebelt groups conducted at room temperature without drying catalyst did not reveal any drying behavior.Results of our own drying experiments with EL at 70 • C, Fig. A3 in SI, are the first published in this field of research that do show drying, as EL conversion goes to completion.The results of our validation are presented in Table 1, showing the set of reaction mechanisms and the associated kinetic rate coefficients.
In this paragraph, we will explain how we obtained the results in Table 1, depending on a) which mechanisms were taken into account and b) for these selected mechanisms, the values of the associated rate coefficients.Note that for all ARNG model simulations with varying numbers of mechanisms and kinetic parameter values, the initial conditions taken were identical to those in Orlova et al. [32], with the exception of the catalysts concentration that we do not include in the system [31].Note that for all ARNG model simulations with varying numbers of mechanisms and kinetic parameter values, the initial conditions taken were identical (except for catalyst, there is none here) to those in Orlova et al. [32]: initial EL concentration 3.71 mol/L, O 2 concentration 2.18⋅10 − 3 mol/L.We choose to scrutinize all the chemical mechanisms that have been proposed in many previous publications.A considerable reduction in the number of rate coefficients was realized by estimating coefficients per reaction family.In previous publications [7][8][9][10]32] distinct values of coefficients were applied for reactions of each species, based on group contribution theory.It turned out that the resulting coefficients are varying by orders of magnitudes, even within the same reaction family.A thorough discussion to explain these differences is beyond the scope of this paper, but this issue certainly needs further investigation.Presently, we decided to refrain from using group contribution theory in an uncritical manner.This implies that we did not use unique coefficients for reactions of each of the thousands distinct monomeric units, but rather have obtained a more general classification of kinetic parameter values for reaction families.Doing so, we were able to construct the relatively concise Table 1, with only 31 coefficient values.Given the limited set of experimental data, it is still expected that while estimating this moderate number of parameters, only few parameters will be found that significantly influence the model outcomes.Generally, the model outcomes that can be experimentally measured are not significantly sensitive to all mechanisms theoretically possible.We identified the 'key' mechanisms for which the model is sensitive, meaning that changing their rate coefficients significantly influence the prediction of experimentally measurable quantities.Note that we have marked those key parameters with an asterisk in Table 1 and that we could estimate them using the ARNG model and experimental data.
In this section we first explain the choice of included mechanisms and then we discuss how this choice influences the number of species produced in the ARNG model.Considering this choice of mechanisms, we propose values or ranges of values of kinetic coefficients.Subsequently, we discuss them in the light of theories relating chemical structure to reactivity, where we take group contribution theory as a reference for comparison.The result of this part are the rate coefficients listed in Table 1, which gives the best possible agreement with the measured FTIR trends of unreacted EL concentration and with SEC-data both monomeric and dimeric speciesof polymerized ML.Furthermore, we performed an additional validation test involving a much larger number of species than the few, based on FTIR and SEC, employed in the parameter estimation procedure leading to Table 1.We compared the ARNG-predicted spectrum to an experimental ESI-MS spectrum of polymerizing EL, properly illustrating the high resolution of our modeling approach, fully matching the fine resolution of mass spectroscopy.To further highlight the results from a polymerization perspective, we present a sensitivity analysis of crosslink types and gel formation to kinetic parameters.The latter property requires the use of RG-theory using ARNG-predicted concentrations as input [31].Finally, we demonstrate the reaction network approach lying at the heart of ARNG, by showing the computed reaction path of one of the acidic groups, as a final product of EL autoxidation polymerization.

Selection of mechanisms
The total number of species and the reactions they participate in was determined in the ARNG model by choosing the reaction rules to be active.We have employed three different sets as shown in Table 2.The question is, are all mechanisms available in the ARNG model really required to provide a good fit with the measured data?Set 1 has the smallest number of reactions and has been constructed to see whether it might generate fair approximative results at low computational expense.Note that it comprises the basic autoxidation mechanisms like Habstraction and OOH-decomposition, in combination with the two β-scission mechanisms.There is a clear difference in model size between Set 1 and Set 2, which is mainly caused by peroxide crosslink decomposition, which is mainly caused by peroxide crosslink decomposition, since implementation of decomposition leads to a larger computation time by a factor of 250.Additional accounting for epoxidation further increased computation time.Fig. 3 shows the results for EL conversion and concentration of C --C bonds of the three sets given the parameter values optimized for Set 2, as compared to FTIR-measurements.Between Sets 1 and 2, the figure shows clear differences, while compared to Set 3, where epoxidation is added, there are no changes in the results.Consequently, we chose to optimize the parameters of Set 2. In this paper we have assumed the absence of a common catalytic drier, such as driers based on cobalt [34].The presence of drier would increase the number of monomeric species even more, but this is beyond the scope of this work.

Kinetic parameter cases for species and reactions
Using Set 2 of reactions to describe the polymerization process, kinetic rate parameters were optimized for validation with experimental data from FTIR spectroscopy.Parameter estimation was done manually using the ARNG model in view of the computationally intensive solution of the model (see Table 2), which does prevent using automatic parameter estimation.This involves manually changing the coefficients to get the best possible fit between the predicted and experimentally observed quantities.We followed an estimation strategy involving a) the set of available experimental data, b) a selection of key reaction mechanisms and c) values of rate coefficients from earlier studies, mainly calculated from structure-reactivity relationships [7][8][9][10][46][47][48]. The quantitative experimental FTIR data used in parameter estimation were the EL concentration, conjugated C --C bond concentration and, to a lesser extent, both the concentrations of carboxylic acid and the sum of alcohol and aldehyde groups.As the FTIR response depends on extinction coefficients of specific groups or species, only the unreacted EL concentration from FTIR is quantitative, up to a certain extent, as the experimental starting concentration is known.For conjugated double bonds the shape of the time profile was used for estimation.Given the set of data, the most influential reactions in predicting measurable quantities are listed in Table 3.The results of these cases for EL, conjugated double bond and hydroperoxyl concentration time profiles are shown in Fig. 4.
The parameters in Table 3 have been determined by adjusting their values such that a good fit is obtained between model and experimental FTIR data for EL concentration and conjugated double bond concentration profile, see Fig. 4.Not unexpectedly, more than one set of parameter values could provide a fair fit with FTIR data.In Table 3 Case I and Case II mark the lower and upper limits of parameter ranges, within which a good fit with FTIR data could be achieved, see Fig. 4. The lower value of hydroperoxyl decomposition in Case II agrees with values predicted by structure-reactivity relations, the higher value in Case I represents the optimum value, given all experimental data, as will be discussed below.Note that many combinations of parameter values within the limiting values of Cases I and II would give a good fit with FTIR-data.However, these values turn out to be highly correlated.For instance, a lower hydroperoxyl coefficient requires higher values for the remaining coefficients for a good fit.Similar strong correlations were not found for other reaction families, like β-scission.
The reaction families listed in Table 3 turned out to be the key mechanisms, which was not surprising in view of following chemical considerations.H-abstraction is important, since hydroperoxides are formed downstream from EL conversion, so increasing the rate of Habstraction does increase the formation of hydroperoxides.Conversely, a slower rate of hydroperoxyl decomposition builds the hydroperoxyl concentration over time due to the mass balance of the formation by peroxyl radical H-abstraction and decomposition into hydroxyl and alkoxyl radicals.Since addition reactions directly occur on conjugated double bonds those groups were expected to be sensitive to addition rate, but their impact on EL concentration is less obvious.
Considering the key processes in Table 3, the values of the kinetic coefficients of H-abstractions are reasonably well in line with previous modeling work [18,[46][47][48][49].They are also in the range of those reported in older sources of kinetic data for peroxyl radical abstraction from a bisallylic H atom [50,51].However, they are 2-3 orders of magnitude lower than what would be expected from structure-reactivity relationships [46][47][48][49]52].For Case I addition rate coefficients are well in line with previous modeling work [18,[46][47][48][49], while also being in the same order of magnitude as values for propagation reactions occurring to single double bonds in radical polymerization.However, for Case II they turn out to be unrealistically high.Finally, hydroperoxyl decomposition in Case I is considerably higher than the band width estimated from earlier work [46][47][48].Fig. 4 shows hydroperoxyl concentration profiles over time, for Case I featuring the maximum that has often been observed in earlier work.Unfortunately, no experimental data for the hydroperoxyl concentration, like peroxide value (PV) [20,53], are available to compare with for uncatalyzed EL drying.Finally, we note that the values of the peroxyl crosslink decomposition rates in Table 3 have been chosen consistently lower than hydroperoxyl decomposition.This is because crosslinks created during drying are mainly peroxyl groups, which are observed to remain stable on much longer time scales [18].
To conclude this paragraph, we note that Table 3 is a significant result, as it demonstrates that possible parameter values of key reaction families leading to a good fit with FTIR data are bound between the limiting Cases I and II.However, it also implies that we need additional information to find the most probable combination of parameter values within these bounds.In the next paragraphs we discuss extra sources of information: acid formation and size exclusion chromatography of methyl linoleate, especially in relation to hydroperoxyl decomposition.

Low hydroperoxyl decomposition rate and presence of carboxylic acid
With Case II we simulated a situation with low values of hydroperoxyl and peroxyl crosslink decomposition.The value of hydroperoxyl decomposition is at the high end of a band width, from 10 − 9 to 10 − 15 s − 1 , based on a study of reaction families relevant for condensed-phase oxidation of hydrocarbons by Pfaendtner and colleagues and a group additivity calculation for the heat of reaction, see Table 4.These low decomposition values are a consequence of the uncatalyzed nature of our system.Cobalt and other metal driers are often used to catalyze the autoxidation process by overcoming the barrier for hydroperoxide decomposition [10].However, for a system without any appreciable drier and the reported inertness of coated rutile TiO 2 [54], the spontaneous fission is not expected to be competitive.Also, other plausible pathways seem to be absent for bimolecular interaction of hydroperoxides with an abstractable hydrogen to yield alkyl and alkoxy radicals and water as products.Hence, structure-reactivity theory apparently points at low hydroperoxyl decomposition.However, in Case II the ARNG predicted amount of carboxylic acid was much smaller than appears from our own FTIR-data, as will be shown below, which casts serious doubt on low hydroperoxyl decomposition.We will discuss the effect of hydroperoxyl decomposition on carboxylic acid and also on hydroperoxyl and alcohol groups as well as ether crosslinks, since some experimental data from SEC [39] and Raman spectroscopy [55] are available that may shed light on this issue.The model also reveals a clear impact of hydroperoxyl decomposition on the alcohol groups and ether crosslinks as Case II has more than two orders of magnitude less OH-groups and O-crosslinks than Case I. Prior Raman data [34] of EL drying at elevated temperature without drier indeed confirm the presence of OH and ether crosslinks.Generally, during the drying of linseed oil or its associated model compounds, the majority of crosslinks produced are peroxyl crosslinks, in addition to a fair amount of ether crosslinks.This is the case for both catalyzed and uncatalyzed processes [34,35,37,56,57].However, even for low hydroperoxyl decomposition detectable amounts of OH-groups and ether crosslinks are formed, so the presence of these groups does not provide us with decisive information about the hydroperoxyl decomposition rate.
This situation is different for carboxylic acid.According to the model, formation of acid may occur through several chemical pathways, one of which is shown in Section 3.5.The model shows that the impact of hydroperoxyl and peroxyl crosslink decomposition on the production of carboxylic acids is much more pronounced than on OH-groups and Ocrosslinks.Fig. 5 shows the concentrations of conjugated double bonds and carboxylic acid groups relative to initial EL concentration on a logarithmic scale.Conjugated double bonds are taken as a reference, since in both FTIR and ARNG model their significant presence is evident.The predicted maximum concentration of carboxylic acid is 12 orders of magnitude lower than the conjugated double bond maximum.In contrast, in the FTIR spectra shown in SI, Fig. A4, carboxylic acid peaks at 1412 cm − 1 and 1715 cm − 1 are clearly visible.Although the extinction coefficient of the carbonyl stretch vibration of carboxylic acids is relatively high (estimated a factor 10 higher than conjugated double bonds at 987 cm − 1 ), this difference can not explain the 12 orders of magnitude lower acid concentration.
We conclude that the presence of carboxylic acids is a strong contraindication of the very low hydroperoxide decomposition in Case II.Also, as remarked before, in order to obtain a good fit with the FTIR data on EL conversion and conjugated double bond concentration, unrealistically high values for addition had to be applied (see Table 3).Consequently, we conclude that the low hydroperoxide decomposition is not in line with available experimental data.We hypothesize that the deviation from structure-recativity theory may be caused by the unintended small amount of a catalyzing substance, like for instance the silica coating on the titanium dioxide [44].

High hydroperoxyl decomposition rate -monomer and dimer mass distribution from SEC
Case I marks the upper limit for the hydroperoxyl decomposition rate.This upper limit has been determined with help of data from size exclusion chromatography (SEC) of a thin film of methyl linoleate (ML) dried at 80 • C for 3.5 h without a catalyst (Mourão [39]).We make use of the fact that hydroperoxyl decomposition leads to the transformation of hydroperoxyl substituted to alcohol substituted compounds having a smaller mass.The ARNG model has been employed to compute the molar mass distribution for monomeric (molecules without a crosslink) and dimeric (two monomeric units connected by one crosslink) species for Case I and for a case with ten times faster hydroperoxyl decomposition, as shown in Fig. 6.The procedure to find the molar masses and concentrations of monomers and dimers is described in the SI, Paragraph A3.
Molar mass distributions shown in Fig. 6 were calculated using the ARNG model for different drying times for the cases discussed in previous paragraphs to arrive at similar relative heights of dominating monomers and dimers: Case I: 16 h and the case with faster hydroperoxyl decomposition: 15 h.Since in SEC concentration peak height is weighted with the square of molar mass, we applied the same weighting to the ARNG-spectrum.The computed mass distribution has a high resolution as the exact molar masses are calculated, while in a SECspectrum peaks are significantly broadened.For clarity and better comparison with SEC artificial peak broadening has been applied to the ARNG mass spectrum, while the total surface under the curve is still normalized to total molar massas in SEC.At the chosen drying times ML conversions are over 95 %, yet unreacted ML is still the most dominant peak.For Case I the main oxidized monomer peak is observed at M = 376 g/mol for monomeric species ML-OH(OOH) 2 and at M = 734 g/mol for the peroxyl crosslinked dimer (OOH) 2 -ML-OO-ML-OH(OOH).When applying a higher hydroperoxyl decomposition the dominant peaks are different: at M = 342 g/mol for monomer ML-(OH) 3 and dimer (OH) 2 -ML-OO-ML-(OH) 2 with M = 686 g/mol.Higher OH-substitution in latter case is attributed to the ten-fold higher hydroperoxyl decomposition rate.This could be confirmed by inspecting the reaction pathway in the model.Hydroperoxyl decomposition leads to O-radicals that perform H-abstraction and this yields OH-groups.
We will compare the ARNG-computed molar mass distribution to the SEC-trace measured by Mourão [39] after 3.5 h drying of ML at 80 o C as shown in Fig. A4 of the SI.Note that the SEC-trace is given with elution time on the horizontal axis, since no calibration to molar mass was done.Mourão attributes the observed peaks to ML, oxidized monomers, dimers and trimers.She interprets the oxidized peaks as hydroperoxyl ML, in line with intuition of oxidized peak identity.However, given the relative distances between the ML peak and the dominant oxidized monomer and dimer peaks, we could confirm that indeed the best match exists with the mostly peroxide-substituted species.To further assess this attribution, calculated higher oligomer masses could have been considered as well, but direct reconstruction of trimers from ARNG as we did for dimers leads to a combinatorial explosion practically impairing such an exercise.In an earlier publication [31] we have applied a coarsegrained method based on Random Graph theory to compute the molar mass distribution for higher oligomers.This approach has allowed us to perform a model-based calibration of the SEC-curve and gives further support to the dominance of peroxy-substitution after the drying times realized in our and Mourão's experiments.On the basis of Fig. 6 and the comparison with SEC we conclude that the Case I value of hydroperoxyl decomposition coefficient, 5⋅10 − 6 L/(mol.s), is indeed the maximum value still complying with the SEC-data.Considering that Mourão's EL samples have been prepared at 80 • C and our kinetic coefficients are estimated at 70 • C, yet a slightly higher maximum value is possible.Assuming an activation energy E a = 40 kJ/mol for hydroperoxyl decomposition, the rate at 80 • C would be a factor 1.5 higher than at 70 • C.However, a full account of the effect of small temperature differences on all reactions is beyond the scope of this paper.Therefore, despite the 10 degrees temperature difference between FTIR and SEC, we maintain our conclusion that the Case I value of hydroperoxyl decomposition coefficient is the best possible estimate.

Conclusion of parameter estimation procedure
In Paragraph 3.2.1 it was concluded that additional information is required to find the most probable combination of parameter values within the bounds marked by Cases I and II.Based on the significant presence of carboxylic acid, we think that the hydroperoxyl decomposition must be considerably higher than predicted from structurereactivity relationships (Case II value).However, given the information from the SEC-trace of ML, it was also observed that the decomposition coefficient should not exceed the Case I value.Therefore, we think that the Case I parameter set has the greatest support from all sources of data, finally leading to the complete set of Table 1.We also have observed that some mechanisms appear to be more 'key' than others.This issue is further illustrated by an extended sensitivity analysis described in Section A4 of the SI, where further model-predicted concentration profiles, for different rate parameter values, compared to experimental (mostly FTIR) profiles of EL, conjugated double bonds, single C --C double bonds and carboxylic acid and oxygen uptake are shown (Figs.A5-A10).

ESI-MS analysis
Molar masses of non-radical, monomeric and dimeric species were obtained from the ARNG model using the same procedure as in Section 3.2.2 and employed to identify, by visual inspection, individual m/z observed in the ESI mass spectrum.Such spectra were collected after varying drying times at 70 • C. To compare ESI with ARNG, we assume that no further polymerization takes place in EL samples in the time interval between drying and subjecting to ESI at room temperature [34].The ESI mass spectrum of EL (both positive and negative modes) has around 820 significant (representing 99.99 % of cumulated peak intensities) peaks at m/z ratios up to 870 g/mol, which is the maximum possible in this model system.ARNG generates 250 monomers (species without crosslinks) up to 406.2 g/mol (species with 3 hydroperoxides, 1 alkene, and an ester) and 56,832 dimers, with the smallest dimer being a species with 144.1 g/mol.
Peaks in the ESI spectrum are identified with the ARNG-predicted monomers.Each peak in the positive or negative ESI mode was associated to a molar mass by subtracting the weight of the attached Na + , Ac − or NH 4 + ion, or adding the mass of a proton to deprotonated molecules.
Thus, each peak has multiple contributions to the mass spectrum.The molar mass of monomers or dimers predicted by ARNG falling in a tolerance of 0.15 m/z around each ESI values was collected.Thus, one ESI product-peak may become matched to one or more model-predicted monomers or dimers.Within the monomer matches, around 300 unique ESI mass peaks out of the total number of 5064 ESI peaks were present in the monomer range, so the largest part of the ESI peaks was attributed to higher oligomers or non-attributable.Likewise, not all ARNG-predicted monomers were found in ESI peaks.The partial overlap of ESI and ARNG also applies to dimers.
The next step in the procedure was assigning the matched peaks ('identified' peaks) to the molar mass of the closest model-predicted molecule.A matched peak may become associated to more than one model species, for instance isomers or near-isomers, very close in molar mass.To solve this problem, the measured peak counts were attributed to model molar mass in proportion to the model-predicted concentrations.As the final step in the procedure the peaks positively identified by the ARNG model as monomers and/or (there may be overlap) dimers were collected in a molar mass histogram with an interval of 0.15 g/mol.Within each interval, the peak intensities, corrected for multiplicity using proportional peak attribution as explained above, were added up.In intervals with both monomers and dimers, the peak count was distributed in proportion to the model-predicted concentrations of those monomers and dimers.
We now conduct a quantitative analysis of the ESI and ARNG-data in the monomer range.Rather than directly comparing all the concentrations either observed in ESI or predicted by ARNG, we will focus on a few monomeric products and checking whether they are co-existing in ESI and ARNG.Fig. 7 shows the ESI-spectrum of relative peak intensities without heating (light grey dots) and after a heating period of 24 h (dark grey dots) and 146 h at 70 • C (black dots).Peaks were normalized with the sum of all monomers peak intensities identified as monomers by the ARNG model, comprising all grey, black and red dots in Fig. 7.Not visible in this plot is the noise level of the normalized ESI peaks at 10 − 4 (see Fig. A11 in SI).There are 42 red dots corresponding to species detected by the ARNG model at non-negligible concentrations, >9⋅10 − 6 mol/L, this is the ESI-ARNG set.More information about these 42 species is available in Table A1 of SI, providing molar masses, ESI peak intensities, ARNG concentrations and the ARNG species identification number.Most of the black dots were present in the ARNG-predicted species, as can be seen from the logarithmic peak intensity plot in Fig. A11 in SI, but within concentrations that we believe too low to be detectable by ESI.In Fig. 7, the 11 numbered red dots are the dominant ESI peaks (>0.05), which were all predicted by the ARNG model.Although not all were dominant peaks.The properties and structural formulas of these 11 monomers are listed in Table A2 of SI.Not all significant ESI peaks are identified as significant by the ARNG model: clusters of non-identified lower ESI peaks appear in ranges 135-165 and 235-265 g/mol, except at 152.1 g/mol.
We consider the ARNG-predicted concentrations and dominance of species in Fig. 8, which has a logarithmic concentration scale to account for the large differences in concentration levels.Again, the red dots mark the monomers that are also found at significant peak intensities in the ESI measurements.The 11 monomers dominant in ESI all appear as significant but not dominant among the ARNG-predicted species.Black dots refer to monomers predicted by the ARNG model, but not detected in the ESI measurements as indistinguishable from noise.For smaller molar masses this might be attributed to the volatility of these species.Most dominant monomers in the ARNG model are also present in the ESI measurement, except for that with the highest molar mass possible: 406.3 g/mol.We have no explanation for this compound not showing up in the ESI measurement.
Further many interesting questions are remaining in the comparison between ARNG results and MS-ESI measurements.One is the possible correlation between the relative magnitudes of the concentrations as predicted by the ARNG model and revealed by the ESI measurement.This is discussed in the SI, Section A5 and Fig. A12.
We made a similar but much more concise analysis of the 56,832 dimeric species.Fig. A13 of SI shows the ESI mass spectrum for the dimer species up to the heaviest dimer at 775 g/mol.Agreement with the model is more complete with the identified species; however, peaks at the low and high end of the dimer range are not identified.This lack of dimer species also suggests greater fragmentation into smaller molecules, such as products deriving from allylic H-abstraction, which is not incorporated in this model.In addition, peaks in the high dimer region (that is dimers without fragmenting) also result from adducts, two monomers "bound" together by the ESI ionizing ion, e.g.NH 4+ ; these are absent in the model.

Crosslink and polymer network formation
Although modeling the polymer network formation is beyond the scope of this work, a brief sensitivity analysis of crosslink and polymer network formation (gelation) to some rate coefficients and a brief discussion of the effect of polymerization on rate coefficients will conclude this sub-section.Polymer network computations have been made by an existing RG-model [31,58].The RG-model uses the connectivity information of all the species generated by the ARNG model as input.Each species may have from zero to three crosslinks, corresponding to 0-to 3functional monomeric nodes, which are mainly peroxyl and ether type crosslinks.From the node-functionality distribution the RG-model computed the global (or polymer) properties, like size distribution and the appearance of an 'infinite' component (a polymer network, or gel).For the purpose of this paper, we only considered gel point and gel fraction.
Fig. A14 of SI shows that under the circumstances simulated (no drying catalyst) peroxyl crosslinks dominate, but still ether crosslinks were present in significant amount, while alkyl crosslinks were absent.This is in line with experimental observations (Oyman et al. [34][35][36]).Fig. A15 of SI demonstrates the strong effect of Russell termination on the distribution of numbers of crosslinks per monomer unit (functionality).Fig. 9 shows the gel or network formation at various peroxyl crosslink decomposition rates for the standard value of Russell termination coefficient (blue, red, green) and for the ten times faster Russell termination (orange).In the latter case no gel is formed at all.
Considering rate coefficients in relation to the progress of polymerization, we note from Fig. A14 (right) that crosslink formation closely follows EL conversion (see Fig. 3), in agreement with development over time of the molar mass distribution found by Van't Hoff et al. [31].Fig. 9 shows the gel point to be at 3 h, while the gel fraction subsequently considerably increases, similar to EL conversion and crosslink formation.Polymer (gel) formation is an important indicator of the liquidsolid transformation during drying.This transformation is expected to cause several 2nd order reactions to become diffusion controlled, resulting in a decrease in those rates.In comparable studies of photopolymerization of multifunctional acrylates we do indeed observe such a decrease [59][60][61][62][63].In the present study rate coefficients were assumed to be constant.From the acrylate studies we speculate that the rate coefficients found in this linoleate study are representative of the early stage of the drying process.Future investigation of drying oils will require to account for the strong negative feedback of polymerization progress on the reaction rates.

Reaction pathway analysis
The ARNG model generates a directed bipartite reaction network of nodes and edges (connections) with 4016 species nodes and 1,552,793 reaction nodes for the Set 2, see Table 2.All the reactions leading to a given species as an (intermediate) product could be identified as the ingoing edges to that species node.Likewise, all reactions consuming a given species were the outgoing edges to all those reactions.This allows us to track and analyze the reaction pathway of any intermediate or end product.As an example, in Fig. 10, we analyzed the reaction pathway of end products hexanoic acid (116.1 g/mol) and 9-ethoxy-9-oxononanoic acid (216.1 g/mol).We also choose these acids to demonstrate the plausibility of the production of acids given the kinetic assumptions, as already discussed in Section 3.2.2.According to the ARNG model these acids are produced in equal amounts (having undergone the same reactive transition, only differing in unreactive tails, SI Fig. A2), each 0.575⋅10 − 4 mol/L in 25 h.Note that these acids are identified in the ESI spectrum, Fig. 7, while the heavier acid features in Table A2 in SI as the rank eight most important monomer.Fig. 10 shows the whole pathway from EL to the acids of interest, through eight intermediate products and nine different reactions.Since the solution of the ARNG model generates the concentration profiles over time of all species, as shown in Fig. A16 of SI, the total amounts of products formed and/or reacted in 25 h could also be computed.These amounts are shown as well in the figure, together with reaction type.The amounts of products are gradually decreasing in downstream direction of the reaction pathway.In Fig. A16 of SI one observes that EL (nr 3), original concentration 3.71 mol/L, fully reacts via H-abstraction (HA) to radical species 5.A large part of species 5 reacts by O 2 addition to hydroperoxyl compound 13, of which 40 % reacts into two steps to the crosslinked group species 284.A small fraction of species 284 undergoes second crosslinking, which undergoes reversible reactions (OOH-decomposition, OH-radical recombination) to species 452.Large part of species 452 (6.89⋅10 − 4 mol/L) undergoes β-scission type 2 yielding aldehyde 531.Further around 20 β-scission type 2 reactions from other compounds also yield 531, in total 1.18⋅10 − 3 mol/L.Just over 1 % of species 53, 1.48⋅10 − 5 mol/L, finally yields the acids through Baeyer-Villiger reactions with peroxyl acids.The most important peroxyl acid is species 1963, responsible for 13 % of the production of hexanoic and 9-ethoxy-9-oxononanoic acid.In this particular pathway, OOH-decomposition is taking place, but peroxyl crosslink decomposition is not.It should be noted, however, that the occurrence of this reaction strongly enhances acid production: by a factor of 7.

Discussion and conclusions
The main achievement of this work is a comprehensive and validated model that describes reactions involved in drying of oil paint.By exploring ethyl linoleate polymerization we realized a major step forward in the kinetic modeling of key steps in the drying of biobased oils, like linseed oil.Our objective was to provide a rigorous quantitative kinetic model using automated reaction network generation, with more predictive power than the numerous qualitative analyses of what happens during drying offered in the past."Oil monomers" are chemically complex as they possess multiple reactive sites of varying propensity on each of its three fatty tails leading to an explosion of possible monomers.Nevertheless, we were successful in validating key the features like EL concentration decay rate and conjugated double bond concentration, featuring a maximum, using FTIR data.We deduced a probable set of kinetic rate parameters (Table 1) for the main autoxidation mechanism in EL polymerization.
Three reaction sets were evaluated with two sets of kinetic parameters estimated from experimental data.A reduced reaction set without epoxide formation (Set 1), for which we could not confirm a discernable effect on, generated over 4000 monomer fragments and 1.5 million reactions between them.Kinetic rate parameters for reaction families such as hydrogen abstraction, O 2 addition, radical termination and addition to conjugated double bonds were manually estimated using the reduced reaction set and FTIR data.However, the range of parameter values that satisfied agreement of the model with species concentrations derived from FTIR spectra had still to be reduced.
We subsequently followed arguments made about acid formation and SEC-data to study the hydroperoxide decomposition rate in more detail.This rate appeared to be much higher than expected in an uncatalyzed system (without drying agent).Our new model revealed that theoretically predicted OOH-decomposition rate would lead to negligible acid concentrations, which were clearly observed in the FTIR spectrum.Also, we could establish a feasible reaction pathway from EL to acids, especially hexanoic acid.We attribute this to unintended catalytic effects of trace components in the oil sample.In addition, we put an upper limit to OOH-decomposition, based on yet another experimental source: a size exclusion chromatography trace of methyl linoleate.Comparing the model outcomes with the SEC-trace and considering the ratio between peroxide-and alcohol-groups allowed us to find an optimal value for OOH-decomposition rate.

Fig. 1 .
Fig. 1.(a) Chemical structure of linoleate (EL or ML, reactive site is highlighted) and its representation as a molecular graph with the abbreviated notation R 1 for the ester part and R 2 for the acid tail (b) The scheme that describes the key steps in the ARNG model (Orlova [32]).

Fig
Fig. 2. Directed and undirected crosslinks on monomeric species.Recombination termination leads to undirected alkyl crosslink edges (blue) on both monomers participating in the reaction.Addition crosslinking leads to an outgoing crosslink edge (red) on the lower molecule and to an ingoing crosslink edge on the upper molecule.Reactive patterns correspond to those found in the Supplemental Information, pattern nr 16 for undirected alkyl crosslinks, nr 43 for ingoing and nr 45 for outgoing OO-crosslinks.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .Fig. 8 .Fig. 9 .
Fig. 7. ESI spectrum of relative peak intensities in the monomer range of EL samples at the start (without heating, light grey dots), after 24 h and after 146 h heating at 70 • C (black).The start is after a few minutes of polymerization, where still EL is detected (308.1 g/mol), while its peak has disappeared from the spectrum at later stages.Red dots mark the ESI peaks predicted by the ARNG model, assuming 146 h of heating at 70 • C, of concentrations larger than 9⋅10 − 6 mol/L.Numbers refer to the height ranking of the 11 most dominant ESI peaks (>0.05), see Table A2, SI.Dominant peaks are all identified by the ARNG model (Table A2, SI).Clusters of lower ESI peaks in ranges 135-165 and 235-265 g/mol are not reproduced by the ARNG model, except at 152.1 g/mol.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .
Fig.10.Analysis of the reaction pathway from EL to hexanoic (116.1 g/mol) and 9-ethoxy-9-oxonanoic acid (216.1 g/mol), total production in 25 h of 1.15⋅10 − 4 mol/L.Numbers refer to the species numbers assigned by the ARNG-algorithm.The T's denote the unreactive end tails and ester tails, as shown in Fig.A2of the SI.The two acids have the same reaction product part (species nr1009), but a different tail.Type of reaction is shown along with total amount reacted from reactant(s) to product(s) in 25 h.Fig.A16of SI shows all the species' concentration profiles over time.EL (species 3), original concentration 3.71 mol/L, fully reacts away by Habstraction(HA) to radical species 5.After a number of reaction steps aldehyde species 531 is produced by β-scission type 2. Just over 1 % of species 531, 1.48⋅10 − 5 mol/L, finally produces the acids through Baeyer-Villiger reactions with all peroxyl acids.The most important acid within this group of peroxyl acids is species 1963, responsible for 13 % of the production of hexanoic and 9-ethoxy-9-oxononanoic acid.

Table 2
Summary comparison of species and reaction sets.

Table 3
Variation of hydroperoxyl and peroxyl crosslink decomposition rate.Numbers refer to reaction families in Table1.