Universal Kinetic Solvent Effects in Acid-catalyzed Reactions of Biomass-derived Oxygenates

SUPPORTING INFORMATION Universal Kinetic Solvent Effects in Acid-catalyzed Reactions of Biomass-derived Oxygenates Theodore W. Walker(a)†, Alex K. Chew(a)†, Huixiang Li(b, c), Benginur Demir(a,d), Prof. Z. Conrad Zhang(b), Prof. George W. Huber(a), Prof. Reid C. Van Lehn(a)*, and Prof. James A. Dumesic(a,d)* (a) Department of Chemical and Biological Engineering, University of Wisconsin – Madison, Madison, WI, 53706 (United States) (b) Dalian National Laboratory for Clean Energy, Dalian Institute of Chemical Physics, Chinese Academy of Sciences, 457 Zhongshan Road, Dalian 116023, Liaoning, China (c) University of Chinese Academy of Sciences, Beijing 100049, China (d) DOE Great Lakes Bioenergy Research Center, University of Wisconsin-Madison, Madison, WI, 53706, United Sates † These authors contributed equally *send correspondence to: vanlehn@wisc.edu and jdumesic@wisc.edu

Section S1: Experimental Details Table S1 presents the kinetic solvent parameters that describe the rates of the seven Brønsted-acidcatalyzed reactions considered in this study, including those data reported in our prior work, 1 in aqueous mixtures of -valerolactone (GVL), 1,4-dioxane (DIO), and tetrahydrofuran (THF).These kinetic descriptors capture changes in the apparent rate constant for each reaction as a function of solvent composition at a fixed temperature, compared to the same rate constant in pure water ( : where ( , ) is the apparent rate constant for the i th reaction, and the subscript denotes the identity and composition (in j th mass fraction) of the organic phase.Herein, we demonstrate that these descriptors are themselves weak functions of temperature and reactor headspace, which allows for direct comparisons to be made between reactions that occur at different temperatures, and in solvents that exert significantly different vapor pressures.Confidence intervals were calculated at the 95% confidence level based on the propagation of error resulting from two known sources of uncertainty, which are discussed in the sections that follow.Details relating to experimental protocols and reaction kinetics modeling are also discussed.
Table S1.Kinetic solvent parameters describing the rates of the Brønsted-acid-catalyzed reactions considered in this study as a function of solvent composition.Reaction conditions: see Table S2.

T rxn = 373 K
Mass fraction (j) Of the organic phase

S1.2 Methods
Reaction kinetics measurements for tert-butanol (TBA), and fructose (FRU) dehydration were carried out using methods described in in prior work. 1 All other reactions were carried out in closed, 10 mL thick-walled glass batch reactors.In a typical experiment, 2 mL of solution, consisting of an appropriate amount of reactant (e.g., xylitol, XYL) and triflic acid catalyst (TfOH), in a solvent system of desired composition (e.g., 90 wt% GVL in water) were filtered using a 0.2 m membrane (VWR International; PTFE) prior to being added into closed batched reactors.
Owing to the limited solubility of some reactants in organic solvents (e.g., XYL), the initial filtration step ensured that only fully solubilized materials were present under reaction conditions, eliminating the possibility of transport effects confounding the kinetic analyses.Concentrations of TfOH in each experiment were varied to maintain consistent kinetic profiles in all solvent systems with respect to reaction time.For example, the concentration of TfOH used for XYL dehydration in pure water was 1.3 M. In contrast, the concentration of TfOH used for the same reaction in 90 wt% GVL/water was 0.037 M. Reactors were placed in an oil bath at the desired temperature and stirred at 500 rpm using magnetic stir bars.The reactors were allowed to equilibrate to the desired reaction temperature for a period of ~15 minutes, after which the reactors were removed at intervals corresponding to the desired reaction times, and quenched in an ice bath at 273 K.
After quenching the reactors, the contents were diluted 9:1 by mass with HPLC grade water, neutralized with sodium bicarbonate to a pH of 7.0, and stirred vigorously for 30 seconds.The neutralized solution was then filtered using a 0.2 m membrane (VWR International; PTFE) prior to analysis.Reaction products for ethyl-tert-butyl ether (ETBE) hydrolysis in water and in GVL/water mixtures were analyzed using a gas chromatograph (Shimadzu GC-2010) equipped with a flame ionization detector.All other analyses were performed using a high-performance liquid chromatograph (Shimadzu LC-20AD) equipped with a differential refractometer (Shimadzu RID-20A) and a photodiode array (Shimadzu SPD-M20A).Concentrations of all analytes were quantified based on calibration curves using external standards.Separation of XYL, 1,4anhydroxylitol, ETBE, ethanol and TBA was achieved using two Bio-Rad Aminex HPX-87P HPLC ion-exclusion columns in series (7.8 x 300 mm, 7m).A mobile phase of HPLC grade water at a flow rate of 0.6 mL min -1 was used.Separation of cellobiose, glucose, and levoglucosan was achieved with an ion-exclusion HPLC column (Bio-Rad; Aminex HPX-87H; 7.8 x 300 mm, 5 m).A mobile phase of 5 mM sulfuric acid aqueous solution (Ricca Chemical; HPLC grade) at a flow rate of 0.6 mL min -1 was used.Separation of ETBE, ethanol and TBA in water in GVL/water mixtures was achieved in a gas chromatography column (Shimadzu; SHRXI-5MS; 30 m, 0.25mm ID, 0.25 mm DF) to establish carbon balances.To minimize exposure of the GC column to mineral acids and other ions, all other reactions of ETBE were monitored using the HPX-87P HPLC system, and the reaction was assumed to be quantitative in these cases as well.

S1.3 Reaction kinetics modeling
Rate constants were estimated in MATLAB (nlinfit function; Levenberg-Marquardt nonlinear least squares algorithm) by optimizing their values to minimize the residuals between experimental and model-predicted time courses in each experiment.Model-predictive expressions in the MATLAB routine took the form of rate laws based on transition state theory, formulated as described below.
We assume that the acid-catalyzed reactions considered in this study occur via a sequence of elementary steps, one of which is rate limiting and passing through a protonated transition state.Under these assumptions, 2 the forward rate of each reaction per unit volume is equal to the concentration of the activated complex (or the transition state) in the rate-limiting step ( , multiplied by the frequency of vibration along the reaction coordinate ( : Transition state theory treats the activated complex as being in equilibrium with the solvated reactants and proton. 3Accordingly, the concentration of the activated complex can be written in terms of the activities of the free proton (aH+), and the reactants and/or products (aR/P) that are consumed or produced in the elementary steps prior to the rate-limiting step: where kB is the Boltzmann constant, T is the temperature, h is Planck's constant, and ‡ is the activity coefficient of the activated complex.K is the product of all equilibrium constants for the elementary steps prior to the rate limiting step, and n denotes the stoichiometry of the species that participate in those elementary steps (being positive for reactants and negative for products).The leading rational term in Equation S3 is the frequency factor that results from standard statistical mechanical treatments of the vibrational mode along the reaction coordinate. 4While the reactant is consumed over the course of the reaction, it remains approximately constant and dilute (~ 0.1 M) throughout the experiments conducted in this study.The activity coefficients of all species represented in Equation S3 may thus be approximated as constants, meaning that the rate expression may be recast as a function of molar concentrations: where , is the apparent rate constant for the i th reaction containing all invariant terms (including activity coefficients) in Equation S3, and the subscript denotes the identity and composition (in j th mass percent) of the organic phase.The rate expression for each reaction in the present study thus takes the form of Equation S4, with as many CR/P terms as there are elementary steps involving the consumption of a reactant, or generation of a product prior to the rate-limiting step.

For acid-catalyzed dehydration reactions
We have shown that the removal of water in the carbenium ion formation step is rate-limiting for reactions of TBA and PDO. 1 It has been shown elsewhere that the first water elimination step is rate-limiting in the acid-catalyzed dehydration of FRU, 5,6 and XYL dehydration involves only a single water removal. 7The rate expressions for these four dehydration reactions are thus: For acid-catalyzed hydrolysis reactions It is generally accepted 8,9 that the addition of water in Brønsted-acid-catalyzed hydrolysis of ethers is fast compared to the formation of the carbocation that follows removal of the leaving group.As such, the concentration of water should not appear in the forward rate expression for Brønsted-acid catalyzed hydrolysis reactions.Accordingly, the forward rate expressions for each of the hydrolysis reactions in this study are: Regarding the reversibility of the LGA hydrolysis reaction In solvent systems containing less than 50 wt% water, glucose undergoes dehydration to afford LGA.
LGA hydrolysis (Equation S11) is thus an equilibrium-limited reaction in such solvent systems, and the reverse rate of this reaction must be accounted for in the kinetic modeling scheme.Accordingly, the forward rate of LGA hydration is measured by accounting for the reaction affinity (or the extent of the system's departure from its equilibrium state) in the rate expression. 10We thus consider the overall rate of the LGA hydrolysis reaction as the difference in the rates of LGA hydrolysis and glucose (GLU) dehydration: where again, the apparent rate constants contain the activity coefficients for the kinetically relevant species, which are approximated as constant.The equilibrium constant for this reversible reaction is: where is the concentration of the i th species at equilibrium.Even in solvent systems consisting of 90 wt% of the organic phase, the concentration of water is in significant excess compared to LGA (~5.6 M compared to ~0.1 M, respectively).As such, is approximately equal to , and combining Equations S12 and S13 affords: Moreover, by definition 11 : (S15) So that, finally, Equation (S11) can be written as: where Korg,j is the equilibrium ratio of the concentrations of glucose and LGA: The value of Korg,j was measured independently of the kinetic experiments in this study by allowing a 0.1 M solution of LGA to equilibrate, in the presence of a triflic acid catalyst, in each of the six solvent systems where the rate of the glucose dehydration reaction competes with the forward rate of LGA hydrolysis.Figure S1 presents the reaction time courses for these six experiments, and captures the relationship between LGA and glucose in terms of the pseudo-equilibrium constant, Korg,j at 403 K. LGA hydrolysis to afford GLU in the presence of a triflic acid catalyst to demonstrate the equilibrium relationship between reactant and product as a function of solvent system.Reaction conditions: 0.3 -0.15 M TfOH; 403 K. Error in the apparent equilibrium ratio of glucose to LGA represents the 95% confidence interval, based on the variability in the two to four data points in each graph after the two species have equilibrated.

Section S2: Experimental Results
Having measured the values of Korg, j, experimental reaction time courses for LGA hydrolysis were fit to Equation S17 in the MATLAB routine, yielding a set of values for the apparent forward rate constant ( . ) in each solvent system.In a similar fashion, Equations S5-S10 were fit to experimental time courses to afford the apparent rate constants for the Brønsted-acid-catalyzed reactions for XYL, CEL and ETBE as a function of solvent composition.Figure S2 contains exemplary reaction time courses for each reaction in pure water, and in 90 wt% GVL, along with the associated model fits and carbon balances.Table S2 presents the apparent rate constants measured in this study, along with the reaction conditions for each experiment.Selectivities differing from unity generally reflect uncertainties in quantification of reactants and products, which is captured in the 95% confidence intervals reported along with each rate constant.In calculating the 95% confidence intervals, we account for the propagation of errors resulting from two sources of uncertainty: uncertainty in the quantification of reactants and products (based on triplicate experiments) and the residuals between experimental and model-predicted time courses (MATLAB; nlparci function).

S2.1 Effects of temperature and reactor headspace on the apparent rate constants
Control experiments were conducted to assess the sensitivity of the kinetic measurements to reaction temperature.Figure S3 demonstrates the weak temperature dependence of the kinetic solvent parameters associated with XYL and PDO dehydration in GVL/water mixtures, indicating that these kinetic descriptors may be compared between reactions that take place at different temperatures.We conducted control experiments to probe the sensitivity of the kinetic measurements to differences in the partitioning of water and the organic cosolvent between the liquid and gas phases.Table S3 demonstrates that even when the most volatile cosolvent (THF) is used, the apparent rate constant for XYL dehydration is essentially invariant with respect to reactor headspace.As such, differences in volatility between the three cosolvents results do not lead to significantly different concentrations of water in the liquid phase under reaction conditions.

Section S3: Molecular Dynamics Simulation Details
Unbiased molecular dynamics (MD) simulations were performed using the CHARMM36 all-atom force field 12 with the TIP3P water model.Molecular parameters not included within CHARMM36 were generated using the CGenFF force field, 13,14 which is fully compatible with CHARMM36.Molecular dynamics was performed using a leapfrog integrator with a 2 femtosecond time step.Verlet lists were generated using a 1.2 nm neighbor list cutoff.Van der Waals interactions were modeled with a Lennard-Jones potential using a 1.2 nm cutoff that was smoothly shifted to zero between 1.0 nm and 1.2 nm.Electrostatic interactions were calculated using the Smooth Particle Mesh Ewald method with a short-range cutoff of 1.2 nm, grid spacing of 0.12 nm, and 4 th order interpolation.Bonds were constrained using the LINCS algorithm.All simulations were carried out in the NPT ensemble using Gromacs 2016 in a cubic simulation cell with periodic boundaries conditions in all directions 15 .All thermostats used a 1.0 ps time constant and all barostats used a 5.0 ps time constant with an isothermal compressibility of 5.0×10 -5 bar -1 (specific thermostats/barostats are detailed in Section S3.1).Simulation configurations were output every 10 ps and the final 190 ns of each production trajectory were used for analysis.

S3.1 System preparation
Each reactant/cosolvent system was prepared using the workflow schematically illustrated in Figure S4.First, water and cosolvent molecules were added to a cubic simulation box with 6 nm box vectors at the desired mass fraction.The method used to calculate the number of cosolvent and water molecules in each system is described in Section S3.2.The solvent mixture was equilibrated for 5 ns at T = 300 K and P = 1 bar with a velocity-rescale thermostat and Berendsen barostat.The reactant was then added to the system and the system was equilibrated for 500 ps at the reaction temperature (see Scheme 1 in main paper) at 1 bar using the same thermostat and barostat.The system was then simulated for 200 ns of production at the same temperature and pressure using the Nose-Hoover thermostat and Parrinello-Rahman barostat.

Figure S4
. Schematic depiction of mixed-solvent system preparation.R denotes the reactant.Note that the second production trajectory, used to calculate hydrogen bonding lifetimes, was not always 4 ns; some reactants required a longer simulation time to obtain accurate hydrogen-bonding lifetime data.

S3.2 Calculation of number of molecules within system
The number of water and cosolvent molecules required for each solvent mixture was determined based on the total volume of the simulation box.The molecular volume of each cosolvent species was estimated by calculating the volume occupied by a single molecule in a cubic box, then multiplying the volume by two (1.5 for water) to ensure the molecules appropriately fit when added to the system.The estimated molecular volumes for water, DIO, GVL, and THF were 0.09 nm 3 , 0.24 nm 3 , 0.30 nm 3 , and 0.16 nm 3 , respectively.The molecular weights for DIO, GVL, and THF were 88.11, 100.067, and 72.11 g/mol, respectively, which were used to convert weight fractions to mole fractions.To find the number of water and cosolvent molecules to add to the system given the mixed cosolvent composition, we performed a volume balance to obtain: where and are the number of molecules of component 1 (water) and 2 (cosolvent), and are molecular volumes of component 1 and 2, is the mole fraction of component 1, and is the total volume.

S3.3 Calculation of accessible hydroxyl fraction ()
The accessible surface area (ASA) of each reactant was calculated using the Shrake and Rupley algorithm 16 as implemented in the MDTraj package.The default 0.14 nm probe radius, which is the van der Waals radius of water, and 960 points were used.The van der Waals radii of all atoms were adjusted according to the values from Bondi. 17The final 190 ns of production data from trajectories containing reactant in pure water were used to compute the hydroxyl fraction (see Figure S4).

S3.4 Calculation of radial distribution function
The reactant-water radial distribution function (RDF) was calculated using the "compute_rdf" function available in the MDTraj package using a bin width of 0.02 nm.The center of mass of each molecule was used to compute intermolecular distances.Periodic boundary conditions were accounted for in all calculations.

S3.5 Estimation of cutoff between local and bulk domains
The cutoff radius between local and bulk solvent domains, r cutoff, was defined for each system as the distance from the reactant where the radial distribution function between the reactant and water molecules reaches unity.These cutoffs were determined by calculating the running average of the radial distribution function starting from r=2 nm, which is a distance large enough for all RDFs to plateau at unity (i.e., bulk behavior).We define r cutoff as the radius for which the difference between the running average and the next point (i.e., ri-1) is greater than 0.015.

S3.6 Robustness of the preferential exclusion coefficient ()
The preferential exclusion coefficient, , is plotted in Figure S5 as a function of the cutoff radius between the local and bulk solvent domains for XYL and TBA in 90% and 50% wt%DIO/water mixtures.For comparison, the cutoff radius for XYL in 90% and 50% DIO/water mixtures is calculated as 1.51 and 0.83 nm, respectively, using the approach in Section S3.5, while the cutoff radius for TBA in 90% and 50% DIO/water mixture is 1.67 nm and 1.37 nm, respectively.The plot illustrates that the  is sensitive to the cutoff radius, but for each system the  plateaus within the error bars by the value of the cutoff radius used in the main manuscript.In Figure S5(b),  for TBA in various 90 wt% organic cosolvents are plotted to show that  converges within the error bars for different cosolvent systems.The cutoff radius for TBA in 90 wt.% in DIO, GVL, and THF is 1.67, 1.57, and 1.73 nm, respectively.Figure S6 presents  as a function of the amount of simulation time used to calculate its value (sampling time) and the initial size of the simulation box.The dependence of  on sampling time was calculated by partitioning the production trajectory in increasing increments of 0.5 ns from the end to the beginning of the trajectory and calculating  for each partition.We find that 95 ns of simulation time is sufficient to determine accurate  values for both TBA in 50 wt% DIO and XYL in 50 wt% GVL.In addition, we find that varying the initial simulation box size does not significantly change the  values.Therefore, we selected 6 nm as the initial box length for all simulations and 95 ns as the simulation block required to calculate reliable  values; note that the box length decreases in each trajectory due to the barostat.

S3.7 Hydrogen bonding lifetimes ()
After the 200-ns production trajectory, an additional 4 ns production trajectory was generated with configurations output every 0.1 ps to obtain hydrogen bonding lifetimes.Hydrogen bonding lifetimes were computed by splitting the 4 ns trajectory into two 2 ns trajectories and calculated using the hydrogen bonding lifetime analysis tool in Gromacs 5.0.1 [18][19][20][21] .Some systems required a longer production trajectory to obtain reliable hydrogen bonding lifetimes based on the standard deviation between the two trajectories.These systems are listed in Table S4.In particular, it was difficult to obtain accurate hydrogen bonding lifetimes for ETBE because it is only a hydrogen bond acceptor.

S3.8 Rescaling of multidescriptor correlation model
Coefficients for Equations 7-9 in the main text were calculated using the "fitlm" function in MATLAB 2017a.The values for Γ, , and were rescaled between 0 and 1 using the following procedure: where is the original simulated descriptor and is the rescaled descriptor.This rescaling procedure allows us to compare the coefficients in the multidescriptor correlation model.

S3.9 Simulation size, cutoff radii for local domain, and data for multidescriptor correlation model
Table S5 presents the number of molecules, simulation box size, and the cutoff radius used to define the local solvent domain for each simulated system.Note that the cutoff radii for pure water are presented as a comparison to the cutoff radii in the cosolvent mixtures, but they were not used to calculate any simulation-derived parameters.

S4.3 Coefficients of multidescriptor correlation model without re-scaling procedure
Coefficients for the multidescriptor correlation model (Equation 9 of the main text) calculated without rescaling molecular descriptors are shown for all cosolvent systems in Table S6.Note that the re-scaling procedure (described in Section S3.8) does not change the accuracy of the multiparameter correlation model.

S4.4 Error in multidescriptor correlation coefficients
The standard error of the best-fit coefficients used in Equation 9 was computed using the "fitlm" function in MATLAB 2017a, shown in Table S7.
Table S7.Standard error in the coefficients for re-scaled and non-rescaled multidescriptor correlation models.
Re-scaled ( where is the variance and RSS is the residual sum of squares.The maximum value is found by differentiating Equation (S21) with respect to the variance, and equating to zero.The result is: The AIC then becomes: For sample sizes where n/k < 40, the AIC is corrected as follows: Corrected AIC (AIC ) values for four different correlation models are shown in Table S8.We find that the four-parameter fit (corresponding to the multidescriptor correlation model in Equation 9of the main text) has the lowest AIC values for DIO and GVL (two out of three cases), indicating that inclusion of and improves the model to a statistically significant extent.

Figure S1 .
Figure S1.Time courses forLGA hydrolysis to afford GLU in the presence of a triflic acid catalyst to demonstrate the equilibrium relationship between reactant and product as a function of solvent system.Reaction conditions: 0.3 -0.15 M TfOH; 403 K. Error in the apparent equilibrium ratio of glucose to LGA represents the 95% confidence interval, based on the variability in the two to four data points in each graph after the two species have equilibrated.

Figure S2 .
Figure S2.Exemplary experimental and model-predicted time courses for reactions of XYL, CB, LGA, and ETBE in water, and in aqueous mixtures of 90 wt% GVL.Reaction conditions: see Table S2.Error bars represent 95% confidence intervals on the data based on triplicate experiments.

Figure S3 .
Figure S3.Effect of temperature on the apparent rate constants for XYL and PDO dehydration in GVL/water mixtures, as expressed by their kinetic solvent parameters ( ).

Figure S5 .
Figure S5. as a function of cutoff radius for (a) XYL and TBA in 90% and 50% DIO/water mixtures.(b)  for TBA in different 90% organic cosolvent mixtures is plotted to show convergence across cosolvents.The dashed vertical lines indicate the cutoff radius used for the calculation of the corresponding  in the main manuscript.

Figure S6 .
Figure S6. as a function of sampling time for (a) TBA in 50 wt% DIO and (b) XYL in 50 wt% GVL for different initial box lengths.The sampling time is presented in increments of 0.5 ns.The values within the parentheses in each legend are the final box lengths after 95 ns of sampling due to the action of the barostat.The dashed line represents the value of  reported in the main manuscript (which is averaged between two independent 95 ns intervals).
-descriptor correlation model for DIO Figure S7 compares values of  calculated using the two-descriptor correlation model (Equation 8 of the main text) to experimentally determined values for DIO/water mixtures.The two descriptors used in this model are the preferential exclusion parameter, Γ, and hydrogen bonding lifetime ratio, .This model is further improved by including the accessible hydroxyl fraction, , as a third descriptor, leading to the multidescriptor correlation model described in Equation 9 of the main text.

Figure S7 .S4. 2
Figure S7.Comparison of kinetic solvent parameters calculated using the two-descriptor correlation model (pred) to experimentally determined values (exp) for seven reactants in DIO/water mixtures.Each reactant has four data points for 0.25, 0.50, 0.75, and 0.90 mass fractions of the organic phase, with the exception of PDO in a 0.75 mass fraction DIO mixture (see TableS1).The slope of the best-fit line for all the data points and the average root-mean-squared error (RMSE) between the values of pred and exp are shown at bottom right.The solid black line indicates a perfect correlation () and dotted lines are drawn at 0 and 0 to help visualize false positive/negative predicted values.Lines above and below the line are shifted by 0.10, denoting the approximate experimental error.

Figure S8 .
Figure S8.Comparison of kinetic solvent parameters calculated using the multidescriptor correlation model (pred) to experimentally determined values (exp) for seven reactants in (a) GVL/water mixtures and (b) THF/water mixtures.Each reactant has four data points for 0.25, 0.50, 0.75, and 0.90 mass fractions of the organic phase, with the exception of reactants in 0.50 mass fraction THF mixtures (see TableS2).The slope of the best-fit line for all the data points and the average root-mean-squared error (RMSE) between the values of pred and exp are shown at bottom right.The solid black line indicates a perfect correlation () and dotted lines are drawn at and to help visualize false positive/negative predicted values.Lines above and below the line are shifted by 0.10, denoting the approximate experimental error.

Table S2 .
Apparent rate constants as a function of solvent system for the four acid-catalyzed reactions examined in this study.Reaction conditions: 2 mL of solution in closed, 10 mL thickwalled glass reactors stirred at 500 rpm.
Accordingly, rate constants for ETBE hydrolysis in these solvent systems were estimated based on the rate of ETBE disappearance only.

Table S4 .
Hydrogen bonding lifetime production trajectories that required more than 4 ns to obtain accurate hydrogen bonding lifetimes.

Table S5 .
Simulation size and number of molecules for the 91 systems studied.Norg and NH2O are the number of organic cosolvents and water molecules, respectively.V is the average volume of the cubic simulation box.rcutoff is the cutoff radius between the local and bulk solvent domains.Γ, , and are descriptors used in Equation 9 of the main text.

Table S6 .
Coefficients of the multidescriptor correlation model (Equation9) calculated without re-scaling molecular descriptors.

Table S8 .
Akaike information criteria with a correction AIC ) for finite sample sizes.