Thermochemical hydrolysis of macroalgae Ulva for biorefinery: Taguchi robust design method

Understanding the impact of all process parameters on the efficiency of biomass hydrolysis and on the final yield of products is critical to biorefinery design. Using Taguchi orthogonal arrays experimental design and Partial Least Square Regression, we investigated the impact of change and the comparative significance of thermochemical process temperature, treatment time, %Acid and %Solid load on carbohydrates release from green macroalgae from Ulva genus, a promising biorefinery feedstock. The average density of hydrolysate was determined using a new microelectromechanical optical resonator mass sensor. In addition, using Flux Balance Analysis techniques, we compared the potential fermentation yields of these hydrolysate products using metabolic models of Escherichia coli, Saccharomyces cerevisiae wild type, Saccharomyces cerevisiae RN1016 with xylose isomerase and Clostridium acetobutylicum. We found that %Acid plays the most significant role and treatment time the least significant role in affecting the monosaccharaides released from Ulva biomass. We also found that within the tested range of parameters, hydrolysis with 121 °C, 30 min 2% Acid, 15% Solids could lead to the highest yields of conversion: 54.134–57.500 gr ethanol kg−1 Ulva dry weight by S. cerevisiae RN1016 with xylose isomerase. Our results support optimized marine algae utilization process design and will enable smart energy harvesting by thermochemical hydrolysis.

Scientific RepoRts | 6:27761 | DOI: 10.1038/srep27761 To convert macroalgal biomass into fuel molecules via a biochemical pathway, the first step is to degrade the cell wall material into fermentable sugars with acid or through enzymatic reactions 18 . However, energy efficient macroalgae cell deconstruction and saccharification into fermentable sugars is still a major challenge 5,15,19,20 . Though acid hydrolysis or pretreatment has been widely used in the degradation process of biomass [21][22][23] , the details of the process parameters and their importance in determining the process output have rarely been reported 24,25 . Understanding of these parameters is important as this will enable the design of energy efficient biorefineries.
To understand how various parameters could influence the cell wall deconstruction efficiency and their comparative significance, we used Taguchi Robust Experiment Design 26 and Partial Least Square Regression. We also used Flux Balance Analysis techniques to model the potential production of ethanol, butanol and acetone from the Ulva biomass hydrolysates.

Materials and Methods
Macroalgae biomass Material. Macroalgae from Ulva genus were collected near the beach of Ramat Aviv, Israel ( Fig. 1), in May 2015. The biomass dried in an oven at 40 °C until constant weight. The dried biomass was made brittle by liquid nitrogen and then grinded into powder manually in a mortar. The Ulva powder was sieved by 30 mesh sieve to make sure all particle sizes are smaller than 0.5 mm. All chemicals and standards were purchased from Sigma-Aldrich (Israel) if not otherwise mentioned. Thermochemical deconstruction. Thermochemical deconstruction was conducted in 10 mL centrifuge tubes (Nalgene ™ Oak Ridge High-Speed PPCO Centrifuge Tubes (Thermo-Fisher Scientific, CA) in autoclave (Tuttnauer 2540MLV, Netherlands ). For each batch, dried samples of Ulva was weighed on analytical balance (Mettler Toledo, Switzerland) Sulfuric acid (Sigma-Aldrich, Israel) was injected into the tube and the mix was vortexed to make the powder well distributed in acid. After deconstruction, the hydrolysates were neutralized by sodium hydroxide (Sigma-Aldrich, Israel). All the solid/liquid ratio, acid concentrations, hydrolysis time and temperature were according to the Taguchi design in Table S1.
Taguchi orthogonal arrays for thermochemical deconstruction. The goal in these series of experiments was to determine the effects of thermochemical hydrolysis process parameters: temperature, treatment time, %Acid and %Solid load on the extraction of monosaccharaides from Ulva biomass. The possible range of process parameters and their combinations is large. Therefore, to decrease the number of experiments, but still be able to evaluate the impact of each parameter independently, we applied the Taguchi's Robust Design method for the experimental design and analysis 26,27 . A key feature of the Taguchi method is to determine the parameters of the controllable process factors with the goal to minimize the impact of uncontrollable factors (noise) in industrial process. The experimental design with orthogonal arrays allows for analysis that prioritizes the comparative impact of the process parameters on the yields.
The experiments, conducted for the L16 orthogonal Taguchi array, which are needed to determine the individual effects of each of the tested parameters on the extraction yields are summarized in Table S1. The hydrolysis (1ml total volume) was conducted in 10 ml tubes (Thermo-Fisher Scientific, CA). The resulting hydrolysates were neutralized by 5 M KOH to pH 7. All experiments were done in duplicate.
In Taguchi design of experiment the best parameter setting is determined using signal-to-noise ratio (R). In our experiments we used "the larger the better" algorithm type. The ratio R is determined independently for each of the process outcomes (OUT) to be optimized. These process outcomes are: concentrations of glucose, rhamnose, xylose, glucuronic acid and total sugars, solution density and %Yield. In the current context, maximizing R corresponds to obtaining the maximum concentration and extraction yields of monosaccharaides. The ratio R of a specific process outcome OUT in experiment j was calculated by: where K is the number of experiments (in our case K = 16); #Reps is the number of experiment repetitions (in our case #Reps = 2) and m Rep is the measurement of the process outcome (OUT) in the specific repetition Rep of experiment j.
Consider a process parameter P (temperature, treatment time, %Acid, %Solid as appears in Table S1). Assume that P has a value of L in n(P,L) experiments (for example, temperature = 100 appears in 5 experiments: P = temperature, L = 100 and n(P,L) = 5). Let J (P,L) be the set of experiments in which process parameter P was applied at level L. Let: ) be the average ratio R for concrete level L of parameter P. The sensitivity (Δ ) of each outcome (OUT) with respect to the change in a parameter P is calculated as:

OUT O UT OUT
Ranking (on the scale of 1-4, where 1 is the highest) was assigned to the process parameters according to the sensitivity ranges obtained.
Partial least squares regression. The ultimate goal of multivariate regression analysis is to create a calibration equation (or set of equations) that, when applied to data of "unknown" samples, measured in the same manner, will accurately predict the quantities of the constituents of interest. The multivariate calibration models were generated using Partial Least Squares (PLS) regression, with the goal of defining a relationship between different process parameters (P i ) and any process output: where OUT is the output (e.g. glucose concentration, rhamnose concentration, xylose concentration, glucoronic acid concentration, total sugars concentration, %Yield or density), A is an empirical coefficient, and P i are the process parameters (temperature, treatment time, %Acid, %Solid, as appears in Table S1). For each output (OUT), a model consisting of four different process parameters (e.g., six models in total) was constructed. First, we investigated the influence of each P on each OUT using the whole set of physical variables. Then, we used Martens' Uncertainty Test (MUT) to identify variables that are most important to predict different sugar concentrations, release yields, and hydrolysate properties. MUT is a significance testing method that can be implemented when cross validation PLS method is applied. It tests the stability of the regression results and the selection of significant P-variables (Unscrumbler, Version 9.1, Camo, Norway) 28 . Cross-validation is the best model and, indeed, the only alternative we have when we lack sufficient samples for a separate test set 28 . Due to the limited number of samples, statistical parameters for the calibration model were calculated by leave-one-out-cross-validation (only one sample at a time is kept out of the calibration and used for prediction). The performance and relevance of PLS models were further evaluated by computing different statistics. The difference between the predicted and measured sugar contents was expressed as a root mean square error of prediction RMSECV (root mean square error of cross-validation). RMSECV is defined as the square root of the average of the squared differences between the predicted and measured values of the validation objects 28 : where OUT m is the parameter of the sample (m), OUT P is the predicted value of the sample and n v is the number of samples in the calibration stage. The predictive capability of all models was compared in terms of the relative standard error for cross-validation sets (denoted as RMSECV (%)).
Additionally, we used the Ratio of Prediction to Deviation (RPD), which is defined as the ratio of the standard deviation of the measured values to the root mean square error of cross-validation (RMSECV) or prediction (RMSEP) 29 . An RPD value below 1.5 was taken to indicate that the model is unusable, between 1.5 and 2.0 indicated that the model has the potential to distinguish between high and low values, and an RPD value between 2.0 and 2.5 was indicative that quantitative prediction is possible. RPD values above 2.5 were considered to indicate that the predictive capability of the model is excellent.
The log10 IR was calculated to normalize the rhamnose, glucose, xylose and glucuronic acid distributions and all further quantitative analyses were made on both, the transformed and untransformed data 28 .Therefore, both, separately for each OUT p , were used as reference values in our PLS modelling.
Carbohydrate composition analysis by High Performance Ion Chromatography. Dionex ICS-5000 platform (Thermo Fischer Scientific, MA, USA) was used to quantify released monosaccharaides in hydrolysate. Carbopac MA1 and its corresponding guard column were used for separation. Electrochemical detector with AgCl as reference electrode was used for detection. A trinary solvent system was used for elution as shown in Table S2. The column temperature was kept at 30 °C and the flow rate was set to 0.4 mL min −1 . Calibration curves were produced for rhamnose, glucose, xylose and glucuronic acid on gradient to determine the concentration of corresponding substances in the hydrolysate. All uronic acid peaks were integrated and calculated accordingly using the calibration curve of glucuronic acid (UA) for estimation.
Total yield was calculate using the following Eq. 7: where m i (μg) is the mass of carbohydrate i in the hydrolysate sample (Total), m dwUlva (mg) is the total weight of the hydrolysed biomass in each sample. The summed carbohydrates were glucose, xylose, rhamnose and glucuronic acid. The concentrations of the rest of the released monosaccharaides were negligible.

Density determination.
To determine the density of hydrolysates we used the resonating micromembranes (RMMs) method. RMMs operate in deposited liquid droplet environment [30][31][32] , characterizing the velocity response in the frequency domain 33,34 . Figure S1 shows a typical response spectrum of a dry RMM device of 500 μ m excited near the first three fundamental modes, denoted by (01), (11) and (21). Downward frequency shifts of specific vibration modes correlate with added droplets mass 30-32 as detailed below.
The resonance frequency in each mode, of indices (mn), is inversely proportional to the square root of effective mass density, which contains contributions from both the solid membrane film and deposited liquid material above it [30][31][32] . The analytical relationship, corresponding to symmetrical geometry of a circular membrane, is: where α (mn) is a geometrical factor associated Bessel function zeros along membrane nodes [30][31][32] , Y = 130 (GPa) is the film Young's modulus, ε = 4•10 −6 is the pre-calibrated tensile stress, and ρ is the effective mass density including both device and deposited droplet contribution. With the known density of Silicon ρ f = 2650 (kg m −3 ) and intrinsic film volume V f = π R 2 h = 0.59 (pL), the deposited liquid is characterized by an average density ρ 1 and a total droplet volume V 1 , defining the composite effective masses density: The downward frequency shift corresponding to an added mass relative to the a dry RMB is estimated using 30 : In differential measurements of frequency shifts corresponding to given modes, the pair of Eqs 9 and 10 enable us to extract both parameters ρ 1 and V 1 in each experiment. Here we chose to work with mode (34), vibrating at f(34) = 82 KHz in the dry settings. The spectral width of each mechanical resonance is proportional both to dissipation and density. Only the average density results are used in the hydrolysates analysis below. This method analyzes internal reaction fluxes based solely on simple physical-chemical constraints, such as reaction stoichiometry and metabolic flux constraints, without requiring exact enzyme kinetic data 35 . This methodology enables the prediction of organism growth rates, as well as the prediction of minimal and maximal production rates of various metabolic compounds, based only on reaction stoichiometry and directionality. FBA-based approaches have a wide range of applications, including phenotype analysis, bioengineering and metabolic model reconstructions [36][37][38][39] .
In our case, we utilized FBA mathematical framework to predict production rates of ethanol, butanol and acetone from the hydrolyzed Ulva biomass. According to Flux Variability Analysis (FVA) approach 40 , we estimate minimal and maximal limits for the flux through the target molecule (ethanol) transporter reaction as follows. First, the maximal possible organism growth rate (MaxGrowth) under given media is estimated in (Eq. 11): (1) 0 where S is the organism-specific stoichiometric matrix (S mr corresponds to stoichiometric coefficient of metabolite m in the reaction r) and  v is a vector of reaction fluxes in the organism (v r is a metabolic flux through reaction r). v Growth is an artificial reaction representing organism growth rate (converting several molecules like amino acids, nucleotides and others into units of biomass) and v media transporters refer to the set of transporters active for the received growth media. Constraint (1) is mass-balance constraint, enforcing the sum of fluxes for all reactions producing each metabolite to be equal to the sum of fluxes consuming it. Constraint (2) (3) is media constraints, limiting the media consumption to 1gDW* h −1 .
In case of MaxGrowth below ε = 10 −5 , the organism is defined as non-viable and it does not produce target molecule. Otherwise, minimal and maximal target molecule (ethanol, butanol, acetone) production rates are estimated in (Eq. 12) under the maximal growth rate constraint (a common assumption in FBA-based simulations): where v target is the flux through the target-molecule transporter.
Specifically, to predict the biofuel molecule production rates, we utilized one-step simulation framework described in previous work 39 with metabolic models of S. cerevisiae 41 (original and modified to incorporate Xylose uptake, simulating RN1016 strain), of E. coli 42 and of C. acetobutylicum 43 . The composition of the Ulva biomass was based on literature 12,13,44 and appears in Table 1. The carbohydrate concentrations were modified for the simulations to those of rhamnose, glucose, xylose and glucoronic acid measured in experiments 1-32 (Table S10). Maximal and minimal ethanol production rates were estimated for 4 organisms: (i) E. coli; (ii) S. cerevisiae WT; (iii) S. cerevisiae RN1016; and (iv) C. acetobutylicum. For C. acetobutylicum, we also investigate maximal and minimal production rates of acetate and butanol.

Results and Discussion
Ulva biomass hydrolysate carbohydrate analysis. First, we quantified the released quantities of major released carbohydrates from the dried Ulva biomass at each experimental condition. The results are shown in Table 2.
We also analyzed the sensitivity of the specific and total sugar release during hydrolysis to the tested process parameters that can be controlled (Figs 2 and 3, Tables 2 and 3). Increase of temperature did not significantly affect the rhamnose, glucose and xylose yields (Fig. 2a,e,i). Shorter time increased rhamnose, glucose and xylose yields (Fig. 2b,f,j). Increasing %Acid to a certain level led to higher rhamnose, glucose and xylose yields, further increase in acid concentration did not increase the yields of these monosaccharide's ( Fig. 2c,g,k). Low and high %Solid led to the increase of the rhamnose, glucose and xylose yields (Fig. 2d,h,l).
The yields of glucoronic acid were not affect by process time and temperature (Fig. 2m,n). As with other hydrolysis products, increase of %Acid increasing %Acid to a certain level led to higher yields. Increasing %Solids increased the yield.
Process temperature and time did not affect the total concentration of released monosaccharaides (Fig. 2q.r). Increase of %Acid increasing %Acid to a curtain level led to higher yields (Fig. 2s). Increasing %Solids increased the yield (Fig. 1t).
Scientific RepoRts | 6:27761 | DOI: 10.1038/srep27761 We also found that changes in temperature did not affect the total yield (as defined in Eq. 7) of carbohydrates (Fig. 3a). Shorter of longer process times led to the increased yields (Fig. 3b). Increase in the %Acid to a certain level increased the total yield, and high %Acid led to reduction of the total yield (Fig. 3c). Low and higher %Solid load led to higher yields (Fig. 3d).
Changes of the process parameters did not significantly affect the hydrolysate density (Fig. 4). The average R values for each process factor (P) appear in Tables S3-S9 (see Supplementary information). The importance of the change in each of the process parameters on Ulva deconstruction appears in Table 4.
The optimum parameters that maximize each of the output factors (OUT) appear in Table 5.

PLS analyses of macroalgae Ulva.
As previously noted, a PLS model constructed for each chemical constituent was first implemented on the whole four physical parameters, and then only the significant parameters were kept and each model was re-assessed. The results are shown in Table 6 (significant variables highlighted by + and bold in the text). Tables 6 and 7 present the results of PLS modelling of different sugar contents. As can be seen, acidity and solidity were found to be important variables for sugar component estimation. Highest PRD values were estimated for glucose, rhamnose and density. RPD value between 1.5 and 2.0 indicate that the model has a strong potential to distinguish between high and low values of Rhamnose and density. For Glucose, an RPD value of 2.02 indicates that a quantitative prediction is possible using acidity and solidity.    Due to the relatively small number of samples and relatively small dynamic range of reference values (e.g. released monosaccharaides concentrations), even extremely small differences between the PLS model and   (Tables 4 and 5), indicating that in the tested range of parameters, %Acid and %Solid play a more important role on the hydrolysis success than the process time and temperature.
Predicted Ethanol production. We predicted possible minimal and maximal ethanol production rates for each of 32 Ulva biomass hydrolysis experiments and for Ulva biomass media without any carbohydrates using  (i) E. coli (Table 8); (ii) S. cerevisiae WT (Table 9); (iii) S. cerevisiae RN1016 (Table 10); and (iv) C. acetobutylicum (Table 11). The carbohydrate content used in the model (gr kg −1 ) in Table 1 was updated based on the results from the current hydrolysis experiments (Table S10) In addition, minimal and maximal possible acetate and butanol production rates were evaluated for C. acetobutylicum (Table 11). For E.coli, the maximum predicted production rate of ethanol is obtained at experiments 8 and 9 (19.459-20.833 gr ethanol kg −1 Ulva) ( Table 8). This corresponds to the highest extraction yields of glucose with thermochemical hydrolysis. For S. cerevisiae WT the maximum predicted ethanol rates are 49.323-51.929 gr ethanol kg −1 Ulva DW and were also observed for simulated fermentation of hydrolysates from experiments 8 and 9 (Table 9). This also corresponds to the highest extraction yields of glucose thermochemical hydrolysis. Previous experimental studies reported on 62 gr ethanol production per kg of Ulva pertusa after enzymatic hydrolyses 45,46 .
From all simulations, the maximum production of ethanol was achieved with Ulva fermentation using S. cerevisiae RN1016 (Table 10). For this organism, in simulations 1, 2, 5, 6, 11, 15 and 16 the predicted ethanol production rates were the same as for media without any monosaccharides (0-8.746 g kg −1 ) derived from hydrolysates. For other simulations we observed an increase in minimal ethanol production rates from 0 up to 54.134 g ethanol kg −1 Ulva (in exp. 8) and an increase in maximal ethanol production rate from 8.746 g ethanol kg −1 Ulva up to 57.500 g ethanol kg −1 Ulva (in exp. 8). The achieved predicted results seem reasonable, since the major carbohydrates source for ethanol is glucose, reaching maximal concentration of 93.53 g ethanol kg −1 Ulva in exp. 8. Also, these findings are supported in previous work 47 , which predicted 160 g of ethanol per kg of Ulva, with glucose concentration of 310 g ethanol kg −1 Ulva (keeping similar glucose-to-ethanol mass ratio of 2:1).
For C. acetobutylicum, in simulations with growth rates not close to zero, the non zero ABE production rate varied from 1.630 to 2.909 gr ABE kg −1 Ulva (Table 11). The ABE yield of predicted ABE fermentation of Ulva    Table 7. Statistical parameters obtained for the cross-validation.
hydrolysate by C. acetobutylicum (calculated by ratio of the sum of ABE products to total sugars extracted with hydrolysis) was 0.02-0.17 range. In the reported experimental ABE, yield of Ulva hydrolysate fermentation by C. acetobutylicum was in the 0.03-0.32 range for hydrolysates with various supplements such as glucose, xylose and nutrients as in CM2 medium 1 . The ABE yield of the hydrolysate without supplements was 0.08 1 . In our simulations, the maximum predicted ABE yield is 23.955-26.830 gr ABE kg −1 Ulva. However, these results were shown for experiments 2, 5 and 15, where the growth rates of the organisms were close to 0 and these results should therefore be treated with care, as they depend on the organism survival at this specific medium. Specifically, in our simulations (Eq. 11) we first maximize organism growth rate, utilizing all available media components in favor of this task and, only then, estimate the molecule production range. This is the reason that media with more sugars does not necessarily lead to higher ABE production but rather to higher organism growth rates. On the other hand, when the amount of monosaccharaides (specifically, glucose) is very low, it acts as a limiting factor to the C. acetobutylicum growth, leading to near-ε (Eqs 11 and 12) growth rates (like in exp. 2, Table 2), the amount of media components not utilized by biomass-constructing reaction (v Growth , Eqs 11 and 12) is relatively high. Therefore, simulations allow molecule production estimation in such low growth rate scenarios. The limitation of this approach is that since ε is arbitrary, configurations with near-ε growth rates may not be viable. Moreover, if minimal molecule (i.e. ethanol) production in these scenarios is zero, ethanol production is not coupled with organism growth (i.e., the organism can grow completely without producing ethanol).   Table 9. Predicted ethanol yields using Ulva hydrolysates fermentation by S. cerevisiae WT, based on BioLego fermentation simulation.