Modelling of the Acetiﬁcation Stage in the Production of Wine Vinegar by Use of Two Serial Bioreactors

: In the scope of a broader study about modelling wine acetification, the use of polynomial black-box models seems to be the best choice. Additionally, the use of two serially arranged bioreactors was expected to result in increased overall acetic acid productivity. This paper describes the experiments needed to obtain enough data for modelling the process and the use of second-order polynomials for this task. A fractional experimental design with central points was used with the ethanol concentrations during loading of the bioreactors, their operation temperatures, the ethanol concentrations at unloading time, and the unloaded volume in the ﬁrst one as factors. Because using two serial reactors imposed some constraints on the operating ranges for the process, an exhaustive combinatorial analysis was used to identify a working combination of such ranges. The obtained models provided highly accurate predictions of the mean overall rate of acetic acid formation, the mean total production of acetic acid of the two-reactor system, and ethanol concentration at the time the second reactor is unloaded. The operational variables associated with the ﬁrst bioreactor were the more strongly inﬂuential to the process, particularly the ethanol concentration at the time the ﬁrst reactor was unloaded, the unloaded volume, and the ethanol concentration when loading.


Introduction
Vinegar production is a biotechnological process essentially involving the biological conversion of ethanol from a given source into acetic acid. Vinegar can be obtained from alcohol, wine, cereals or fruits, among other sources [1][2][3][4]. The key step of the process is possibly that by which a complex microbiota of acetic acid bacteria (AAB) convert ethanol into acetic acid in a bioreactor. The bacterial mixture affecting the conversion arises from the natural microbial selection in the acetification medium. In practice, only AAB can exist in an environment containing medium concentrations of ethanol and acetic acid at the beginning but low levels of the former and high levels of the latter at the end [5]. These conditions make it unnecessary to sterilize containers or keep aseptic conditions during operation.
Acetification bioreactors usually operate in a semi-continuous mode. Thus, once the reactors are in full operation, the ethanol concentration is allowed to decrease to a preset level and then an also preset fraction of the reactor contents is unloaded, the remainder being allowed to stand in it in order to act as an inoculum in the next conversion cycle [6][7][8][9][10]. After the bioreactor is unloaded, it is slowly replenished with a fresh alcoholic substrate to start a new ethanol depletion cycle.

1.
Once the ethanol concentration in the first reactor was reduced to E u1 , a volume fraction of culture medium V u1 was unloaded, whereas the remainder was used as inoculum for the fresh medium to be loaded in order to replace the unloaded fraction. Immediately before the first reactor was partially unloaded, the second was completely unloaded in order to receive the portion withdrawn from the first.

2.
Both reactors were loaded with fresh substrate (wine) controlling the ethanol concentrations to E l1 in the first and E l2 in the second. This allowed drastic changes in the ethanol and acetic acid concentrations to be avoided. In fact, intermittently adding appropriate amounts of fresh substrate allowed a constant concentration of ethanol to be maintained throughout the loading stage. 3.
Once loading was finished, both reactors continued to operate until the second was again unloaded prior to the first.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 24 involved on the acetification process has been used, considering the fulfilment of several restrictions arising from the operation of the two bioreactors, which required identifying and examining the impact of such variables (specifically, establishing their lower and upper limits in the framework of the experimental design). Once the experiments were conducted, polynomial models of relevant target variables of the process such as productivity, acetification rate, etc. were fitted using the gathered experimental data. As far as we know, this thorough study where a detailed analysis considering the constraints for the operational variables, the so many experiments, as well as the replications carried out for each experimental set of variables' values has not been done before; additionally, the modelling approach for this acetification set-up has not been reported in any other previous work.

Operating Mode
A diagram showing the operation of two bioreactors working serially is depicted in Figure 1, where the process involved the following operations: 1. Once the ethanol concentration in the first reactor was reduced to 1 , a volume fraction of culture medium 1 was unloaded, whereas the remainder was used as inoculum for the fresh medium to be loaded in order to replace the unloaded fraction. Immediately before the first reactor was partially unloaded, the second was completely unloaded in order to receive the portion withdrawn from the first. 2. Both reactors were loaded with fresh substrate (wine) controlling the ethanol concentrations to  This scheme allowed the following operational variables to be changed: ethanol concentration at the time the first reactor was unloaded (E u1 ), volume unloaded from the first reactor (V u1 ), the ethanol concentration while the first and second reactor were loaded (E l1 and E l2 , respectively), and the temperature in each reactor (T 1 and T 2 , respectively). A total of six operational variables were thus involved. Although the alcoholic strength of the wine and the airflow rate could also have been used as variables, the feed wine usually contains 12% (v/v) ethanol, and the airflow rate is usually as high as possible (especially when productivity is to be maximized). Also, using gas condensers minimized the sweeping of volatiles.
Except for the temperatures, the working range for the previous variables could not be freely chosen because each variable was subjected to physical limitations, especially because the serial operation of the two reactors imposed constraints arising from the mutual relationships between the variables. This required careful planning of the experimental work in order to define the number of experiments needed to obtain accurate information with a view to predicting the influence of each variable on the resulting system behavior. As described in Section 3.1, a fractional factorial design for this purpose was used [36]. Figure 2 shows the two serially arranged bioreactors used to conduct the experiments. The reactors were two 8 L Acetators from Heinrich Frings GmbH & Co. KG (Rheinbach, Germany) equipped with self-aspirating turbines in the bottom in order to obtain a dispersion of fine air bubbles. The reactors were highly efficient in transferring oxygen between the gas phase and the substrate [11]. The ensemble was interfaced to a computer for automated operation, especially with regards to loading and unloading of reactors, and monitoring of the process. This resulted in highly reproducible loading and unloading, and in efficient acquisition of data [6,11,37]. This scheme allowed the following operational variables to be changed: ethanol concentration at the time the first reactor was unloaded ( 1 ), volume unloaded from the first reactor ( 1 ), the ethanol concentration while the first and second reactor were loaded ( 1 and 2 , respectively), and the temperature in each reactor ( 1 and 2 , respectively). A total of six operational variables were thus involved. Although the alcoholic strength of the wine and the airflow rate could also have been used as variables, the feed wine usually contains 12% (v/v) ethanol, and the airflow rate is usually as high as possible (especially when productivity is to be maximized). Also, using gas condensers minimized the sweeping of volatiles.

Experimental Set-Up
Except for the temperatures, the working range for the previous variables could not be freely chosen because each variable was subjected to physical limitations, especially because the serial operation of the two reactors imposed constraints arising from the mutual relationships between the variables. This required careful planning of the experimental work in order to define the number of experiments needed to obtain accurate information with a view to predicting the influence of each variable on the resulting system behavior. As described in Section 3.1, a fractional factorial design for this purpose was used [36]. Figure 2 shows the two serially arranged bioreactors used to conduct the experiments. The reactors were two 8 L Acetators from Heinrich Frings GmbH & Co. KG (Rheinbach, Germany) equipped with self-aspirating turbines in the bottom in order to obtain a dispersion of fine air bubbles. The reactors were highly efficient in transferring oxygen between the gas phase and the substrate [11]. The ensemble was interfaced to a computer for automated operation, especially with regards to loading and unloading of reactors, and monitoring of the process. This resulted in highly reproducible loading and unloading, and in efficient acquisition of data [6,11,37].

Microorganisms
The inoculum used was a mixed culture of acetic acid bacteria (AAB) obtained from bioreactors previously operating in repeated acetification cycles with the same type of wine as the substrate. This is the usual mode of operation for the vinegar industry.

Analytical Methods
The only variable not determined automatically was acidity, which was measured by acid-base titration with an NaOH solution approximately 0.5 N that was previously standardized with potassium hydrogen phthalate. The volume and ethanol concentration were measured in a continuous manner by using an EJA 110 differential pressure probe from Yokogawa Electric Corp. (Tokyo, Japan) and an Alkosens probe equipped with an Acetomat transducer from Heinrich Frings, respectively.

Polynomial Models and Fitting Methods
The models used here were based on second-order non-linear polynomial equations of the following type: where Y denotes the dependent variable, b 0 the independent regression term, b i the coefficients of the linear (first-order) terms, b ij those of non-linear (interaction) terms, b ii those of quadratic terms, X i the independent variables, X j the independent variables with i < j, and n the number of independent variables. The Y variable is the response or output of the polynomial model and X variables are the inputs. The equations were fitted by using one of several specific methods along with experimental data to calculate the previous coefficients. Such methods allow identifying those terms in Equation (1) that are significant and those that are not. In practice, they solve a least-squares problem by minimizing the sum of squares of the residuals (prediction errors). Second-order polynomial models can be fitted by using several methods [38][39][40][41][42][43]. In this work, the Forward Stepwise Regression method has been used, which successively incorporates one by one those independent variables that contribute to predicting a dependent variable from the most to the least predictive. The sequence in which terms are added to or removed from a polynomial in each step of the process is established according to the F-to-enter, F-to-remove (reference values set by the modeler for the F-values associated to each term of the model), R and R 2 criteria. Then, the F-statistics also give a measure of the model sensitivity to each term of the polynomial relating the input variables.

ANalysis Of VAriance (ANOVA)
Before the model equations were fitted, the experimental data to be used for this purpose were checked for statistical differences by analysis of variance (ANOVA), which compares the means for two or more data samples in terms of their variances. ANOVA's null hypothesis is that all means are identical and the alternative hypothesis that at least one mean will be different from all others [36]. ANOVA is used to identify sources of differences among data samples and to assess whether the differences among sample media are too large to be assigned to random errors alone [44].

Experimental Design
The experimental design that was initially intended to use was a Box-Behnken or Doehlert central composite factorial design. With two levels per factor, the total number of experiments needed was 2 n + 2n + 1, where the first term represents the overall factorial design, the second the central points of the faces, and the third the central point. The factors (operational variables) used were the ethanol concentration at the loading stage in the first bioreactor (E l1 ), the ethanol concentration at unloading time in the first bioreactor (E u1 ), the volume unloaded from the first bioreactor into the second one Appl. Sci. 2020, 10, 9064 6 of 23 (V u1 ), the working temperature in the first bioreactor (T 1 ), the ethanol concentration at loading time in the second bioreactor (E l2 ), and the working temperature in the second bioreactor (T 2 ). A total of 77 different experiments were in theory needed for to model these operational variables. In practice, however, this would have involved too many runs-probably more than needed to fit the model. Also, each experiment would have to be replicated many times to obtain reproducible results and lag phases would be needed each time the operational variables were modified to investigate a new case.
It was therefore necessary to reduce the number of experiments without sacrificing experimental information. This was accomplished by using a fractional factorial design requiring only 2 n−p experiments, where n is the total number of variables and p the number of those variables obtained from some interactions of the others. Thus, considering p = 2 (variables T 1 and T 2 ) and design generators I 1 = V u1 ·E u1 ·E l1 and I 2 = V u1 ·E u1 ·E l2 for a 1/4 fraction design, the experimental setup is shown in Table 1. Normalized values −1 and +1 have been used for V u1 , E u1 , E l1 , E l2 and columns for T 1 , T 2 have been obtained using I 1 , I 2 generators, respectively. For example, normalized values for T 1 were the product of normalized values of V u1 , E u1 and E l1 ; similarly, normalized values for T 2 were the product of normalized values of V u1 , E u1 and E l2 . Table 1. Normalized values of the operational variables for the proposed 2 6−2 fractional factorial experimental design.
It was thought essential to use more than one central experiment (level 0 of all operational variables) at different times in order to examine the variance of the response variables. In this work, we used two central experiments (see Table 2). We would also use far-spaced experiments (levels −2 and +2) symmetrically placed at a distance α from center, α being the fourth root of the number of fractional factorial design experiments excluding central points (i.e., α = 4 √ 16 = 2). Therefore, the distance of levels (−2) and (+2) from level (0) was twice that from (0) to (−1) and (+1). This required adding two experiments per operational variable (see Table 3). In summary, a total of 30 experiments would be required with the proposed design (i.e., less than one-half those required by a pure factorial design), with the added advantage that the ranges of operating conditions would be extended to far extreme values. Table 2. Normalized values of the operational variables for the central points of the proposed 2 6−2 experimental design. Table 3. Normalized values of the operational variables for the extended points of the proposed 2 6−2 experimental design.
Because of the serial operation of the two bioreactors, several specific constraints must be considered for previous variables, since they cannot take certain values in practice: Constraints on E l1 : 1.
The ethanol concentration of the feed wine would invariably be 12% (v/v).

2.
E l1 ≥ E u1 as it would make no sense to have an ethanol concentration at the reactor loading stage lower than that at unloading time.

3.
E l1 must be less than the maximum ethanol concentration in the first bioreactor (E l1 ) max at any time. The latter concentration depends on E u1 , V u1 and the ethanol concentration of the feed wine, and can be calculated from the following mass balance: Constraints on E u1 : 1. E u1 must be greater than 1% (v/v) in order to prevent ethanol depletion from happening too early in the first bioreactor-the second bioreactor was used for this purpose. Also, too low a value would slow down acetification and decrease productivity as a result [15].

2.
E u1 must not be too high. Otherwise, ethanol not consumed in the first bioreactor would have to be used in the second and the probability of complete depletion before unloading would be diminished. Also, ethanol concentrations higher than 5-6% (v/v) could detract from microbial activity [15,45]. Therefore, E u1 must not exceed 5% (2)).
Constraints on V u1 : This variable must range from 1 to 7 L if the self-aspirating turbine with which the bioreactors were equipped is to operate properly. The turbine was used to feed oxygen to the microorganisms and help homogenize the culture medium.
Constraints on T 1 and T 2 : The values of these variables must be compatible with the activity of acetic acid bacteria (AAB), which are the microorganisms affecting the process. The optimum temperature for AAB to oxidize ethanol into acetic acid is 25-35 • C. Because AAB activity drops outside this range [46,47], temperatures from 26 to 34 • C have been used here.
The ethanol concentration in the feed wine is typically 12% (v/v). 2. E l2 must be less than the maximum ethanol concentration in the medium of the second bioreactor at any time. Such a concentration depends on E u1 , V u1 and the ethanol concentration of feed wine, and can be easily obtained from a simple mass balance: Based on the previous constraints, E u1 , V u1 , E l1 and E l2 are mutually related; also, feasible E u1 and V u1 values span the range 1-5% (v/v) and 1-7 L, respectively. However, not all potential combinations of the values of the operational variables can be used in practice. Thus, the initial volume of the second bioreactor-unloaded from the first-V u1 , together with E u1 , imposes some additional constraints on E l1 and E l2 . As a result, the experimental design is subjected to the following constraints: • Based on constraint 2 on E l1 and the experiments of Table 1, E l1 level (−1) must not be less than E u1 level (+1). Likewise, E l1 level (−2) must be greater than or at least equal to E u1 level (0) (see Table 3).

•
Based on Tables 2 and 3, E l1 level (0) can only be used in conjunction with E u1 levels (0), (−2) and (+2). Therefore, E l1 level (0) must be greater than or equal to E u1 level (+2), unless no experiments under the extended conditions are to be conducted-in which case E l1 level (0) must be greater than or equal to E u1 level (0).
Based on the foregoing, the constraints on E l1 and E l2 can be summarized as follows: Constraints on E l1 : • Level (+1) must be less than the smallest (E l1 ) max value (Equation (2)) when V u1 and E u1 levels (−1) and (+1) are to be used; also, level (−1) must be greater than or equal to E u1 level (+1). • Level (0) must be less than the smallest (E l1 ) max value (Equation (2)) when V u1 and E u1 levels (0), (−2) and (+2) are to be used, but greater than or equal to E u1 level (+2). • Level (+2) must be less than (E l1 ) max (Equation (2) when V u1 and E u1 are used at level (0) and level (−2) must be greater than or equal to E u1 level (0).
A systematic analysis of all possible combinations of the values of the operational variables with provision for the above-described constraints was done by using a self-developed script in MATLAB [48] (see MATLAB file "Get_feasible_combinations.m" in Supplementary Materials). Each combination corresponded to a complete experimental design, as depicted in Tables 1-3. The software used different values for level (0) of each variable within a preset range in combination with different values for levels (−2) and (+2)-levels (−1) and (+1) were calculated automatically in each case. The specific ranges considered for level (0) were 1-5% (v/v) for E u1 , E l1 and E l2 ; and 1-7 L for V u1 , all in 0.5 steps. A 0.5 step over the range 0.5-5 was used for iteration of levels (−2) and (+2) in each case.
The script initially identified 148611 feasible combinations. Such a large number required additional constraints to be imposed based on the following practical arguments: A1. There must be a simultaneous maximum difference between levels (−1) and (+1) in E u1 and V u1 (see Table 1). This constraint expanded the feasible range of operational conditions and allowed the influence of each variable and mutual relations to be more accurately detected and examined as a result.
A2. Because in the semi-continuous operation mode it is a common practice to unload at least half of the working volume from the first bioreactor, V u1 level (−1) must be greater than or equal to 4 L.
A3. Because the initial ethanol concentration in the second bioreactor must not be too high in order to facilitate appropriate performance (see constraint 2 on E u1 ), E u1 level (+2) must be less than 5% (v/v).
Introducing these additional constraints in the MATLAB script reduced the number of feasible combinations to 310 (see tab "Pass 1" in Excel file "Results.xlsx" in Supplementary Materials). Subsequently, in line with the previous argument A1, the MATLAB script was used with those combinations maximizing the difference between E l1 levels (−2) and (+2), which further reduced the number of combinations to 31 (see tab "Pass 2" in Excel file "Results.xlsx" in Supplementary Materials). Finally, selecting the combinations resulting in the greatest differences between E l2 levels (−2) and (+2) reduced the number of choices to 4 (see Table 4). On the other hand, too low an E l2 value would slow down the process exceedingly at the loading stage. This fact allows the combinations involving the lowest E l2 value at level (−2) (viz., 1, 3 and 4 in Table 4) to be discarded, thus leaving a single one which fulfils all constraints imposed (viz., combination 2 in Table 4). Tables 5-7 show the set of experiments corresponding to such a combination as established according to the 2 6−2 experimental design previously described by using an identical range for both working temperatures (T 1 and T 2 ): 26-34 • C.
As recommended when using experimental design methodology, only the experiments in Tables 5  and 6 would initially be performed (in random order). If the results were not conclusive enough, then the experiments in Table 7 would also be needed. Table 5. Experiments performed according to the proposed 2 6−2 fractional factorial design.  Table 6. Experiments for the central points of the proposed 2 6−2 fractional factorial design.  Table 7. Extended experiments for the proposed 2 6−2 fractional factorial design.

Experimental Results
Conducting the 18 experiments of Tables 5 and 6 (compiled in Table S2.1 in file "S2.docx" in Supplementary Materials) required changing the operating conditions between experiments, so cycles had to be repeated until steady conditions were reached with each new set of operation variables; after that, the experiments were replicated to obtain statistically valid results. Table 8 shows the total number of cycles and that of useful cycles for each experiment. As can be seen, the minimum number of cycles used to calculate the target variables was 5-some experiments required nearly 20, however. The fact that the total number of cycles differed markedly among experiments was a result of being performed in random order. Therefore, depending on the differences in the operating conditions between consecutive experiments, either a lower or a greater number of cycles may be required for the system to become stable again as a result. By way of example, Figure 3 shows the results in terms of ethanol, acidity, volume, and total and viable cells in the first reactor in Experiment 1. Each result was the mean of 5 replicate cycles with its standard deviation (see Table 8). Similarly, Figure 4 shows the ethanol, acidity, and volume results for the second reactor as the means for 5 replicated runs and their standard deviations. The previous results were used to obtain the values of Table 9. The procedure used to calculate the values of non-measurable variables (e.g., mean rate of acetic acid formation, total acetic acid production and mean volume) is described in file "S1.docx" in Supplementary Materials. Figure S1.1 in such file shows the typical time course of ethanol concentration, acidity and volume in the first bioreactor and Figures S1.2-S1.3 show the time course of such variables in the second bioreactor when a fast loading or pre-depletion step exists, respectively. Also, the results obtained in the 18 experiments of Tables 5 and 6 Tables 10 and 11 show the experimental values of the key variables intended to be modelled using polynomial models relating them to the operational variables.

Obtained Polynomial Models
Below are stated the polynomial models (or response surfaces) corresponding to the target variables. By way of example, the specific procedure used to estimate the mean overall rate of acetic acid formation in the proposed two-bioreactor system, (r A ) global est , is detailed below, and that for each of the other variables in file "S3.docx" in the Supplementary Materials.
3.3.1. Mean Overall Rate of Acetic Acid Formation ANOVA exposed the presence of significant differences at a 99.9% confidence level in the mean overall rate of acetic acid formation, (r A ) global , among experiments. In the absence of such differences, (r A ) global would be independent of the operational variables due to experimental error. Therefore, the 18 experiments of Tables 5 and 6 were representative enough of the operating conditions, so the other 12 potentially required (experiments in Table 7) were unnecessary. As shown below, the ANOVA led to identical conclusions for all target variables. Accordingly, the experimental (r A ) global data in Table 10 were fitted to a second-order polynomial by Forward Stepwise Regression, as described in Section 2.6.1, using the software SigmaStat 2.0 [49]. Although 27 terms (viz., the 6 individual operational variables and their 21 possible interactions) were initially considered, the maximum possible number was that of the experiments performed: 18.
Each fitted polynomial was checked for significant differences in the regressions between each step and the previous one by the effect of the addition or removal of terms. ANOVA revealed the presence of significant differences at a 95% confidence level in all cases and, therefore, the results are worth no additional detailed description here.
The fitting steps using the Forward Stepwise Regression method (Section 2.6.1) are detailed in Section S3.1 in file "S3.docx" in Supplementary Materials (Tables S3.1 to S3.14 show intermediate results). The final model was that of Equation (4) and allowed the mean overall rate of acetic acid formation to be estimated with an error of 0.01 g acetic acid·(100 mL·h) −1 . Such a model, however, only held over the operating ranges shown in Tables 5 and 6 for each variable.
Based on statistical significance, not all operational variables in Equation (4) had a direct influence on (r A ) global (in fact, only E u1 had a direct impact) and this variable was independent of T 2 . Also, only some of the interaction terms influenced (r A ) global .

Variable Value
Total cycle duration (h) 39.7 ± 0. 0.11 ± 0.01 Mean overall rate of acetic acid formation in the two bioreactors as a whole (% w/v·h −1 ) 0.13 ± 0.01 Total acetic acid production in the first bioreactor (g acetic acid·h −1 ) 10.1 ± 0.3 Total acetic acid production in the second bioreactor (g acetic acid·h −1 ) 6.6 ± 0.5 Total acetic acid production in the two bioreactors as a whole (g acetic acid·h −1 ) 16.7 ± 0.5 Table 10. Experimental values of selected variables: (r A ) global is the mean overall rate of acetic acid formation in the two-bioreactor system, P m the total acetic acid production in the two-bioreactor system, E u2 the ethanol concentration at the time the second bioreactor was unloaded, V u2 the volume unloaded from the second reactor, t cycle the duration of a cycle, and V m the mean overall volume in the two-bioreactor system during a cycle.

Exp. No.
(r A ) global , g Acetic Acid (100 mL·h) −1 P m , g Acetic  Figure 5, which is the same as Figure S3.1 in in file "S3.docx" in Supplementary Materials) reveals that the two coincided exactly in 10 of the 18 experiments and differed only very slightly in the other 8. Also, all predictions fell within the 95% prediction interval. influence on ( ) (in fact, only 1 had a direct impact) and this variable was independent of T2. Also, only some of the interaction terms influenced ( ) . A comparison of the ( ) predictions with the experimental values ( Figure 5, which is the same as Figure S3.1 in in file "S3.docx" in Supplementary Material) reveals that the two coincided exactly in 10 of the 18 experiments and differed only very slightly in the other 8. Also, all predictions fell within the 95% prediction interval.

Figure 5.
Plot of (r A ) global est against (r A ) global , and curves for the 95% confidence and prediction intervals.
As can be seen from Figure 6 (same as Figure S3 As can be seen from Figure 6 (same as Figure S3.2 in in file "S3.docx" in Supplementary Material), the residuals (viz., the differences between experimental and predicted values) had a zero or near-zero mean in most experiments and were normally distributed in all. The process used with the other estimated variables is described in file "S3.docx" in Supplementary Material. Their polynomial models are discussed below.

Total Acetic Acid Production
This variable, Pm, was also estimated from the experimental data of Table 10. ANOVA revealed statistically significant differences at a 99.9% confidence level between experiments and, hence, Pm was dependent on the operational variables. Therefore, the experimental data of Pm were fitted to a second-order polynomial by using Forward Stepwise Regression to construct the model of Equation (5), which estimated this variable with an error of 0.7 g acetic acid·h -1 . The fitting steps are detailed in Section S3.2 in file "S3.docx" in Supplementary Material (Tables S3.15 to S3.28 show intermediate results) Figure 6. Residuals of the fitting of (r A ) global est for each experiment.
The process used with the other estimated variables is described in file "S3.docx" in Supplementary Materials. Their polynomial models are discussed below.

Total Acetic Acid Production
This variable, P m , was also estimated from the experimental data of Table 10. ANOVA revealed statistically significant differences at a 99.9% confidence level between experiments and, hence, P m was dependent on the operational variables. Therefore, the experimental data of P m were fitted to a second-order polynomial by using Forward Stepwise Regression to construct the model of Equation (5), which estimated this variable with an error of 0.7 g acetic acid·h −1 . The fitting steps are detailed in Section S3.2 in file "S3.docx" in Supplementary Materials (Tables S3.15 to S3.28 show intermediate results) As can be seen, P m est depends directly on V u1 , E u1 and E l1 , as well as on various interaction terms. A comparison of P m est and its experimental counterpart ( Figure S3.3 in file "S3.docx" in Supplementary Materials) revealed the presence of scarcely significant differences. Also, the residuals of the fitting had a near-zero mean in most experiments and were normally distributed in all (see Figure S3.4 in file "S3.docx" in Supplementary Materials).

Final Ethanol Concentration at the Time the Second Reactor Was Unloaded
As with the previous two variables, the ANOVA on E u2 (Table 10) revealed the presence of statistically significant differences at a 99.9% confidence level between the predicted and experimental results. As before, the experimental results were fitted to a second-order polynomial by using Forward Stepwise Regression. The fitting steps are detailed in Section S3.3 in file "S3.docx" in Supplementary Materials (Tables S3.29 to S3.41 show intermediate results). The resulting model (Equation (6)) predicted the ethanol concentration with an error of 0.3% (v/v).
As can be seen, the variable E u2 est is directly dependent on E u1 and on some interaction terms; however, it is independent of T 1 .
As with the previous variables, a comparison of E u2 est and its experimental counterpart ( Figure S3.5 in file "S3.docx" in Supplementary Materials) confirmed that the predictions of the model were quite accurate; also, they fell within the 95% prediction interval except in one case.
Finally, as can be seen from Figure S3.6 in file "S3.docx" in Supplementary Materials, the residuals had a near-zero mean in most experiments and were normally distributed in all.

Volume of Fermentation Medium Unloaded from the Second Reactor
The ANOVA of the volume unloaded from the second reactor, V u2 , revealed the presence of statistically significant differences between experiment means at a 99.9% confidence level (see Table 10).
Like the previous variables, V u2 was fitted by Forward Stepwise Regression. The fitting steps are detailed in Section S3.4 in file "S3.docx" in Supplementary Materials (Tables S3.42 to S3.53 show intermediate results). The resulting model (Equation (7)) predicted it with an error of 0.22 L.
As can be seen, V u2 est is independent of E l1 and T 1 . A comparison of V u2 est and its experimental counterpart, V u2 , confirmed the goodness of the predictions (see Figure S3.7 in file "S3.docx" in Supplementary Materials). In fact, all fell within the 95% prediction interval. Also, the residuals had a zero or near-zero mean in most experiments and were normally distributed in all (see Figure S3.8 in file "S3.docx" in Supplementary Materials).

Total Cycle Duration
As in the previous cases, the ANOVA on t cycle revealed the presence of statistically significant differences at a 99.9% confidence level (see Table 10). Applying Forward Stepwise Regression to the experimental data provided the model of Equation (8), which predicted t cycle with an error of 2.5 h. The fitting steps are detailed in Section S3.5 in file "S3.docx" in Supplementary Materials (Tables S3.54 to S3.60 show intermediate results) As can be seen, t cycle est is influenced by the variables V u1 , E l1 , E u1 and T 1 , which is logical since the overall behaviour of the first bioreactor dictates when the second is to be unloaded. A comparison of t cycle est and its experimental counterpart, t cycle , revealed that all predictions fell within the 95% interval (see Figure S3.9 in file "S3.docx" in Supplementary Materials). Finally, as with the previous variables, the residuals had a near-zero mean in most cases and were normally distributed in all (see Figure S3.10 in file "S3.docx" in Supplementary Materials).

Mean Overall Volume in the Two-Bioreactor System
The ANOVA on the experimental data of V m (Table 10) revealed statistically significant differences between experiment means at a 99.9% confidence level. The multiple regression model obtained (Equation (9)) estimated the mean overall volume in each cycle with an error of 0.33 L. The fitting steps are detailed in Section S3.6 in file "S3.docx" in Supplementary Materials (Tables S3.61 to S3.74 show  intermediate results). Also, as can be seen from the equation, V m est is independent of T 1 .
As shown by a plot of V m est against the experimental data V m ( Figure S3.11 in file "S3.docx" in Supplementary Materials), all predictions fell within the 95% interval, so the degree of fitting was acceptable. This was further confirmed by the residuals for each experiment ( Figure S3.12 in file "S3.docx" in Supplementary Materials), which had a near-zero mean in virtually all experiments and were normally distributed.

Mean Ethanol Concentration in the First Bioreactor
The ANOVA of EtOH m1 (Table 11) also exposed statistically significant differences between experiment means at a 99.9% confidence level. The experimental data were fitted by Forward Stepwise Regression to construct Equation (10), which reproduced them with an estimation error of 0.2% (v/v). The fitting steps are detailed in Section S3.7 in file "S3.docx" in Supplementary Materials (Tables S3.75 to S3.81 show intermediate results) As expected, EtOH m1 est is independent of T 1 , E l2 and T 2 . Comparing EtOH m1 est and its experimental counterpart EtOH m1 ( Figure S3.13 in in file "S3.docx" Supplementary Materials) revealed that all predictions of the model fell within the 95% interval. Therefore, the fitting was good, as further confirmed by the residuals (Figure S3.14 in file "S3.docx" in Supplementary Materials), which had a near-zero mean in most experiments and were normally distributed.

Mean Ethanol Concentration in the Second Bioreactor
Based on the results of the ANOVA on EtOH m2 (Table 11), there were statistically significant differences between experiment means at a 99.9% confidence level. Fitting the experimental results by Forward Stepwise Regression provided the model of Equation (11), which estimated EtOH m2 with an error of 0.4% (v/v). The fitting steps are detailed in Section S3.8 in file "S3.docx" in Supplementary Materials (Tables S3.82 to S3.88 show intermediate results) As can be seen, EtOH m2 est is independent of E l1 and T 1 . A plot of EtOH m2 est against experimental data EtOH m2 ( Figure S3.15 in file "S3.docx" in Supplementary Materials) revealed that the former all fell within the 95% prediction interval. Also, the residuals ( Figure S3.16 in file "S3.docx" in Supplementary Materials) were all normally distributed and had a near-zero mean in most cases.

Mean Acetic Acid Concentration in the First Bioreactor
Based on the results of the ANOVA on HAc m1 (Table 11), there were statistically significant differences between experiments at a 99.9% confidence level. Fitting the data provided the model represented by Equation (12), which predicted the mean acetic acid concentration in the first bioreactor with an error of 0.2% (w/v). The fitting steps are detailed in Section S3.9 in file "S3.docx" in Supplementary Materials (Tables S3.89 to S3.95 show intermediate results) As can be seen, HAc m1 est depends on the same variables as EtOH m1 est and is also independent of T 1 , E l2 and T 2 -a logical result, since the two variables are mutually related.
A comparison of HAc m1 est and its experimental counterpart HAc m1 ( Figure S3.17 in file "S3.docx" in Supplementary Materials) revealed that the predictions of the model all fell within the 95% interval. Also, as can be seen from Figure S3.18 in file "S3.docx" in the Supplementary Materials, the residuals had a near-zero mean in virtually all cases and were normally distributed.

Mean Acetic Acid Concentration in the Second Bioreactor
As with the previous variables, the ANOVA on HAc m2 (Table 11) exposed statistically significant differences at a 99.9% confidence level in mean acetic acid concentration between experiments. Equation (13) Similarly to HAc m1 est and EtOH m1 est , HAc m2 est is dependent on the same operational variables as HAc m2 . A comparison of HAc m2 est predictions and experimental data (HAc m2 , Figure S3.19 in file "S3.docx" in Supplementary Materials) revealed that the former invariably fell within the 95% interval. Also, the residuals ( Figure S3.20 in file "S3.docx" in Supplementary Materials) had a near-zero mean in most cases and were normally distributed in all. Therefore, the model can be deemed acceptably accurate.

Discussion about the Obtained Polynomial Models
Although a black-box model does not allow one to ascertain why some operational variables are influential whereas others are not, it could be interesting to identify the most influential variables and their interaction terms. The influence (statistical significance) of each term in a polynomial equation can be assessed through statistic F, which was used here to decide whether a term was to be included or excluded. By way of example, Table S3.14 in file "S3.docx" in Supplementary Materials reveals that the highest F values were those for E u1 and E 2 u1 . Therefore, the variable (r A ) global est was especially sensitive to the ethanol concentration at the time the first bioreactor was unloaded-it was directly influenced by E u1 and by its quadratic term.
As a rule, the operational variables associated to the first bioreactor were more markedly influential on most of the dependent variables than were those pertaining to the second. This is unsurprising if one considers that the first reactor not only contributed to the total acetic acid production but also supplied the second with the microorganisms which must work under more extreme conditions in the second bioreactor, since one of the main goals was to deplete ethanol in the medium. Therefore, the conditions prevailing in the first reactor should allow a high concentration of very active acetic acid bacteria to be maintained. As stated in the introduction, such conditions are obtained by keeping the ethanol and acetic acid concentration at not too high levels. It is thus unsurprising that the polynomials used to estimate (r A ) global est and P m est were so strongly dependent on E u1 and E l1 as the latter two variables are directly related to acidity in the reaction medium. V u1 is also highly influential; in fact, the greater the volume unloaded into the second bioreactor is, the more marked will be the potential changes in ethanol and acidity levels in the first as a result of the need for a greater volume of fresh medium for replenishment. As expected, the interaction term E u1 ·E l1 in the polynomial for EtOH m1 est is especially important; in fact, changes in E u1 and E l1 must have a strong impact on the mean ethanol concentration in each transformation cycle in the first bioreactor.
One other interesting inference is that the polynomials for EtOH m1 est and HAc m1 est are complementary; in fact, they only differ in the independent term and in the signs of the others. The sum of the independent term coincides with the overall content of the medium [% (v/v) ethanol + % (w/v) acetic acid]. Since the total concentration remains constant, in the absence of volatile losses by sweeping-which was the case with our experiments-this result is unsurprising and provides support for the correlation procedure used to develop the equations. Similar reasoning can be applied to EtOH m2 est and HAc m2 est .
The variable E d2 est is of special interest as it is a measure of ethanol depletion in each biotransformation cycle. One aim of the acetification process may be not to operate at the highest possible rate but rather to deplete or nearly deplete the substrate in each cycle-in which case E d2 est will be zero or near-zero. Again, the variables E u1 , E l1 and V u1 will be especially influential-not directly, but through their interaction terms-, but so will T 2 and E l2 .
Once the previously described models have been obtained, the operating conditions can only be optimized, for specific purposes, through a well-designed optimization process using several objective functions. Hence, additional work would be necessary in this regard.

Conclusions
Despite the broad available experience and technical knowledge available on the biotransformation of ethanol into acetic acid in the vinegar production process, a number of essential questions remain unanswered. Such is the case, for example, with the nature of the microbiota that effects the process and with its complex metabolic interactions. In practice, the process continues to require more or less extensive modelling in order to relate operational variables to specific objective functions.
Black-box models based on generalized polynomials have proved especially suitable for representing the behaviour of acetification systems in the form of response surfaces.
In this work, an effective experimental design based on useful data for the unequivocal calculation of the coefficients of the polynomial equations has been developed. Once the main operational variables have been identified, their admissible ranges have been established. That information has been used into an experimental design in order to identify the combinations of values of the operational variables that would allow the number of experiments maximizing the predictive ability of a model to be minimized. Another aim carried out in this work was to model key variables related to the acetification process performed with two serial bioreactors working in the semi-continuous mode by using quadratic polynomial equations.
The operational variables considered were the ethanol concentration at the time of unloading (E u1 ), unloaded volume (V u1 ), ethanol concentration during loading (E l1 ) and operating temperature (T 1 ) in the first reactor, and the ethanol concentration during loading (E l2 ) and operating temperature (T 2 ) in the second. As it has been shown, the variation ranges for these variables are subject not only to physico-chemical constraints, but also to others arising from the fact that the two reactors operate in a serial mode and from the limited number of combinations allowed for by the experimental design.
As a result, a fractional factorial design with 30 experiments (see Tables 5-7) considering the previous constraints has been obtained, with it being necessary to perform only 18 of them (see Tables 5 and 6).
With the gathered experimental data, second-order polynomials for several target variables of interest in terms of the considered operational variables involved in the industrial production process were developed. Specifically, models for the mean overall rate of acetic acid formation in the two-bioreactor system [(r A ) global ], total acetic acid production in the system (P m ), ethanol concentration at the time the second bioreactor was unloaded (E u2 ), the volume unloaded from the second bioreactor (V u2 ), cycle duration (t cycle ), the mean overall volume in the two bioreactors during a cycle (V m ), the mean ethanol concentration in the first bioreactor during a cycle (EtOH m1 ), the mean ethanol concentration in the second bioreactor during a cycle (EtOH m2 ), the mean acetic acid concentration in the first reactor during a cycle (HAc m1 ) and the mean acetic acid concentration in the second reactor during a cycle (HAc m2 ) were obtained. The experimental results and their estimates were correlated via Forward Stepwise Regression, which allowed models with high predictive ability and minimal errors to be obtained. The resulting goodness of fit allowed a set of polynomials to be established that accurately reproduced the experimental results with only 18 experiments rather than the 30 needed in theory. Such models should allow us to carry out further optimization studies.
With all models, the variables associated to the first bioreactor were the more influential on the process. This was particularly so with E u1 , E l1 and V u1 , and can be ascribed not only to the fact that these variables depend on the fermentation conditions in the first bioreactor-and hence contribute to the total production of acetic acid-but also to the operating conditions in the second-usually more extreme-being influenced by those under which the first is operated.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-3417/10/24/9064/s1. "Get_feasible_combinations.m": MATLAB script for systematic analysis of feasible combinations of operational variables. File "Results.xlsx": Excel file with successively obtained feasible combinations of operational variables. File "S1.pdf": Description of the procedure used to determine the non-measurable variables of the process. File "S2.pdf": Description of the experimental results. File "S3.pdf": Detailed description of the procedure used to obtain the polynomial models of all analysed variables.