Optimization 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 wine acetiﬁcation, previous works concluded that using a single bioreactor hindered simultaneously reaching high productivities with high substrate consumption and the use of two serially arranged bioreactors (TSAB) could achieve such goal. Then, the aim of this work is the optimization, using Karush–Kuhn–Tucker (KKT) conditions, of this TSAB using polynomial models previously obtained. The ranges for the operational variables leading to either maximum and minimum mean rate of acetiﬁcation of 0.11 ≤ ( r A ) global ≤ 0.27 g acetic acid · (100 mL · h) − 1 and acetic acid production of 14.7 ≤ P m ≤ 36.6 g acetic acid · h − 1 were identiﬁed; the results show that simultaneously maximizing ( r A ) global and P m is not possible so, depending on the speciﬁc objective, different operational ranges must be used. Additionally, it is possible to reach a productivity close to the maximum one (34.6 ≤ P m ≤ 35.5 g acetic acid · h − 1 ) with an almost complete substrate use [0.2% ≤ E u 2 ≤ 1.5% ( v / v )]. Finally, comparing the performance of the bioreactors operating in series and in parallel revealed that the former choice resulted in greater production.


Introduction
Bioprocesses, in general, and acetification as an example, are especially complex owing to the large number of variables involved and their mutual interactions. When a bioprocess is microbially driven-which is the case with acetification-most of the complexity arises from the biological activity of a cell population that is usually highly diverse and capable of adjusting to a variety of environments-and hence of behaving in rather different manners. The problem is even greater when the microbial population comprises not a single species but rather a mixture of species or even genera. Although omics techniques have considerably increased the available knowledge on industrial bioprocesses, there is still a long way to go before they can be elucidated at the molecular scale, so their modelling is an invaluable tool to design them as efficiently as possible.
Because such systems are normally industrially or economically significant, developed models are frequently used to identify the operating conditions that will optimize the output in terms of specific performance indicators or variables. The scientific literature abounds with very recent references to modelled bioprocesses: the xanthan biosynthesis by Xanthomonas campestris ATCC 13951 using winery wastewater [1], the uricase production by Aspergillus welwitschiae strain 1-4 [2], the use of polyethylene terephthalate (PET) for the cultivation of Pseudomonas sp. GO16 [3], the extracellular protease production by a native Bacillus aryabhattai Ab15-ES [4] and many other examples could be found. Usually, however, the large number of variables involved in a process, and their mutual interactions, require careful analysis and specific optimization methods with a view to establishing their optimum values for the intended purpose. 2 of 16 The acetification process, which involves a bio-oxidative transformation effected by a mixture of acetic acid bacteria (AAB) [5], leads to the biological conversion of ethanol into acetic acid. Industrial acetification bioreactors usually operate in an automated repeated semi-continuous mode, so that once the fermenters are fully loaded, the ethanol concentration is allowed to decrease to a preset level and then a certain fraction of the reactor volume is unloaded, actuating the remainder as an inoculum for the next conversion cycle [6][7][8][9][10][11]. After the bioreactor is unloaded, it is slowly replenished with fresh alcoholic substrate to start a new ethanol cycle. For any specific substrate, typical operational variables are the ethanol concentration at the time the reactor is unloaded, the percentage of unloaded volume and the loading rate [12][13][14][15][16][17], which influence the mean ethanol and acetic acid concentrations in the culture medium leading to more or less stressing environmental conditions for the acetic acid bacteria (AAB) and so affecting to the rate and efficiency of the process. In practice, it is necessary to optimize the outcome of the acetification process under operationally restricted conditions, for instance, the end product must have a very little ethanol concentration but if ethanol in the culture medium is almost depleted, the high acidity and lack of substrate would affect very negatively the action of AAB in next cycles. In this regard, the use of two identical serially arranged bioreactors working in a repeated semi-continuous mode, see Figure 1, could be a proper alternative to overcome the problem; this possibility is also suggested from several reported models for the acetification process [7][8][9][10][11][18][19][20][21][22][23][24][25][26][27]. In many vinegar plants the main reactors (bioreactor B1 in Figure 1) are partially unloaded in additional depletion fermenters (bioreactor B2 in Figure 1) where the ethanol concentration could be exhausted. This set up, in practice, results in six operational variables the values of which must be found: E u1 , ethanol concentrations at the time the first bioreactor is unloaded; E l1 and E l2 the ethanol concentrations during loading of the first and second bioreactor, respectively; V u1 volume of medium unloaded from the first bioreactor and T 1 and T 2 the constant working temperatures of the first and second bioreactor, respectively. The acetification process, which involves a bio-oxidative transformation effected by a mixture of acetic acid bacteria (AAB) [5], leads to the biological conversion of ethanol into acetic acid. Industrial acetification bioreactors usually operate in an automated repeated semi-continuous mode, so that once the fermenters are fully loaded, the ethanol concentration is allowed to decrease to a preset level and then a certain fraction of the reactor volume is unloaded, actuating the remainder as an inoculum for the next conversion cycle [6][7][8][9][10][11]. After the bioreactor is unloaded, it is slowly replenished with fresh alcoholic substrate to start a new ethanol cycle. For any specific substrate, typical operational variables are the ethanol concentration at the time the reactor is unloaded, the percentage of unloaded volume and the loading rate [12][13][14][15][16][17], which influence the mean ethanol and acetic acid concentrations in the culture medium leading to more or less stressing environmental conditions for the acetic acid bacteria (AAB) and so affecting to the rate and efficiency of the process. In practice, it is necessary to optimize the outcome of the acetification process under operationally restricted conditions, for instance, the end product must have a very little ethanol concentration but if ethanol in the culture medium is almost depleted, the high acidity and lack of substrate would affect very negatively the action of AAB in next cycles. In this regard, the use of two identical serially arranged bioreactors working in a repeated semi-continuous mode, see Figure 1, could be a proper alternative to overcome the problem; this possibility is also suggested from several reported models for the acetification process [7][8][9][10][11][18][19][20][21][22][23][24][25][26][27]. In many vinegar plants the main reactors (bioreactor B1 in Figure 1) are partially unloaded in additional depletion fermenters (bioreactor B2 in Figure 1) where the ethanol concentration could be exhausted. This set up, in practice, results in six operational variables the values of which must be found: Because of the previous comments and with the aim of optimizing the process, a first work was carried out for modelling the behavior of two serially arranged bioreactors working in a repeated semi-continuous mode [27]. Now, in this work, the obtained black box models are being used to find the operating conditions leading to the maximum and Among the different approaches used for modelling this process, the black-box models based on second order generalized polynomials [7,11,[24][25][26][27] showed to be the best alternatives because of their simplicity and accurate predictions.
Because of the previous comments and with the aim of optimizing the process, a first work was carried out for modelling the behavior of two serially arranged bioreactors working in a repeated semi-continuous mode [27]. Now, in this work, the obtained black box models are being used to find the operating conditions leading to the maximum and minimum values of specific industrially useful variables that were used as objective functions.

Raw Material and Microorganisms
As substrate, white wine from the Montilla-Moriles D.O. (Córdoba, Spain) containing (11.5 ± 0.5) % (v/v) ethanol and an initial acidity of (0.4 ± 0.1)% (w/v) as acetic acid was used. Imitating industrial procedures, the inoculum used was a natural mixed culture [5] maintained and stored in our laboratories from experiments in a long-term fully operational bench acetator working with either wine or a synthetic ethanol medium; the original inoculum was taken, approximately a couple of years before, from an industrial tank in full operation (UNICO Vinagres y Salsas, S.L.L., Córdoba, Spain). First, a stage for reactivating and for the adaptation of the inoculum, which includes several previous acetification cycles to achieve repetitive results, was performed [27]. Additional details about some methodological aspects for obtaining the experimental data used to build the models necessary for the optimization study carried out in this work, were described elsewhere [27]; in any case, as a summary, 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.

Operating Mode
Operation of the two bioreactors working serially, see Figure 2, is summarized as follows: once the ethanol concentration in the first reactor decreased to a preset value E u1 , a volume of medium V u1 was unloaded into the second. Then, both bioreactors were loaded with fresh wine to an also preset level (E l1 and E l2 , respectively). Once the maximum volume of medium with which each bioreactor could be loaded (8 L) was reached, the acetification system entered an ethanol depletion stage. When the ethanol concentration again fell to E u1 , the second bioreactor was completely unloaded for replenishment with medium from the first. Depending on the particular operating conditions, the second bioreactor may not be loaded to full capacity and V u2 be less than 8 L as a result. In addition, each reactor can be operated at a different temperature (T 1 and T 2 , respectively). Different values of each of the previous variables can therefore lead to minimum or maximum levels of the dependent variables.
One should bear in mind that the primary aim was to maintain the first bioreactor under non-stressing conditions for the culture to retain its activity; in this way, unloading a fraction of fermentation medium into the second bioreactor would supply it with healthy bacterial biota capable of operating under the substrate depletion conditions they would inevitably find in that reactor. By way of example, the profiles of some state variables of the system, namely, volume of medium, ethanol concentration and acetic acid concentration in each bioreactor are shown in Figure 2

Optimization Method
The values of the operational variables providing the optimum polynomial functions were determined by solving the constrained non-linear optimization problem defined by Equation (1): where ( ) denotes the objective function, the vector of operational (decision) variables, ℎ( ) the set of equality constraints and ( ) that of inequality constraints.

Optimization Method
The values of the operational variables providing the optimum polynomial functions were determined by solving the constrained non-linear optimization problem defined by Equation (1): where f (x) denotes the objective function, x the vector of operational (decision) variables, h(x) the set of equality constraints and g(x) that of inequality constraints.
If these functions are assumed to be differentiable, then optimal points fulfilling the Karush-Kuhn-Tucker (KKT) conditions can be determined [28,29] from the Lagrange function shown in Equation (2) for a maximization problem (for a minimization problem, the Lagrange function is the same as Equation (2) simply by changing the sign of f (x)).
where λ j and µ i are the KKT multipliers. Equation (3) sets the first-order KKT conditions needed for optimality: where ∇ x is the gradient operator with respect to x. ∇ x g i and ∇ x h j must be linearly independent of the active constraints at the critical points x* obtained as solutions to Equation (3). The active constraints at x* are those with associated non-zero KKT multipliers. In a maximization problem, if f (x) is concave and all g(x) functions are convex, then x* will be the global solution to the optimization problem. Otherwise, it will be necessary to check whether the second order KKT conditions for x* are fulfilled in order to confirm that they provide a local maximum on f (x). A function is concave if its Hessian matrix is negative semi-definite (i.e., if all its eigenvalues are equal to or less than 0) and convex if it is positive semi-definite.
Sufficient second-order KKT conditions can be verified by examining the sign of the last n − m leading principal minors of the bordered Hessian B (Equation (4)) evaluated on a critical point.
where n is the number of operational variables, m that of active constraints at the critical point, ∇ 2 xx L(x, λ, µ) the Hessian matrix of L(x, λ, µ) and g(x) the set of active constraints. The kth leading principal minor of B(x, λ, µ) equals B(x, λ, µ) k×k [i.e., the determinant of the k × k submatrix taken from the upper-left of B(x, λ, µ)]. Therefore, the last leading principal minor of B(x, λ, µ) equals B(x, λ, µ) (n+m)×(n+m) [i.e., the determinant of the whole matrix B], the penultimate one being B(x, λ, µ) (n+m−1)×(n+m−1) and so on.
For a critical point to be a local maximum, sign B(x, λ, µ) (n+m)×(n+m) must be equal to (−1) n and the previous n − m − 1 leading principal minors must have an alternating sign.

Optimization of the Mean Rate of Acetic Acid Formation in the Two-Bioreactor System
One potentially useful variable for assessing the "health" of a bioprocess and its productivity is the rate of biotransformation. Because the bioreactors operated in a nonsteady state here, the target variable was the mean reaction rate, which can be estimated from the rate of ethanol consumption or that of acetic acid formation [6]. When volatile losses are negligible, as it is the case [27], both rates are identical. We focused on the mean rate of acetic acid formation, (r A ) global , because it is easier to estimate.
In [27], the polynomial model for the mean rate of acetic acid formation in the twobioreactor system was obtained (Equation (5)). The model allowed (r A ) global est to be estimated with an error of 0.01 g acetic acid·(100 mL·h) −1 .
As Equation (5) shows, not all the operational variables have a direct influence on (r A ) global est ; additionally, as explained in [27], particularly in its Supplementary Material ("S3.docx" file), not all the terms in Equation (5) have the same influence on (r A ) global est ; the higher the value of the statistic F (see Table 1), the more significant the term is for (r A ) global est ; it can be seen from Table 1 that the most influential term is that of E u1 .  [27], Supplementary Material ("S3.docx" file)).

Terms
F Coefficient This model is of especial practical interest as it allows the operating conditions leading to maximum-desirable-or minimum levels of the industrially useful variables to be identified. We used the procedure described in Section 2.3 to develop a MATLAB script [30] that allowed the maximum or minimum value of (r A ) global est to be calculated (see file "Optimal_rA.m" in Supplementary Materials).
The optimization problem addressed to maximize (r A ) global est is represented by Equations (6) and (7). The latter sets the constraints imposed by the experimental conditions on the operational variables.
where V u1 is the volume of medium unloaded from the first bioreactor into the second, E l1 the ethanol concentration during loading of the first bioreactor, E u1 the ethanol concentration at the time the first bioreactor is unloaded, T 1 the temperature in the first bioreactor and E l2 the ethanol concentration in the second bioreactor during loading. Solving the previous optimization problem provided the operating conditions of Table 2.  Because of the number of involved operational variables, it is not straightforward to graphically show the optimum obtained with the conditions of Table 2. Nevertheless, considering the previous comments on the significance of the terms in Equation (5), when plotting the response surfaces of (r A ) global est varying just a pair of the operational variables at once, it is possible to show that significant changes can only be found when the variable E u1 is included. By way of example, response surfaces varying E u1 and E l2 on one side and E l1 and E l2 on the other are shown in Figure 3a,b, respectively. Values from Table 2  plotting the response surfaces of ( ) varying just a pair of the operational variables at once, it is possible to show that significant changes can only be found when the variable 1 is included. By way of example, response surfaces varying 1 and 2 on one side and 1 and 2 on the other are shown in Figure 3a,b, respectively. Values from Table 2 are used for the remaining variables.  Table 3. With a 0.2 interval for each variable, such ranges led to a total of 162 combinations of which only 109 (see Table S1 in file "S1.docx" in Supplementary Materials) fell within the previous range of ( ) . However, using all 162 combinations ensured that this variable would be very close to its maximum value: 0.25-0.27 g acetic acid·(100 mL·h) −1 . Therefore, any of the 162 combinations would, in theory, lead to the highest possible mean rate of acetic acid formation in the two-bioreactor system -and hence to values of each operational variable falling in the ranges of Table 3.  Table 4 (see file "Opti-mal_rA.m" in Supplementary Materials). The problem was like that stated in Equations (6) and (7) except that the aim was to minimize ( ) rather than maximize it.   . Response surfaces of (r A ) global est varying E u1 and E l2 (a) as well as E l1 and E l2 (b), using values from Table 2 for the remaining operational variables.
Since (r A ) global est could be subject to an error of up to 0.01 g acetic acid·(100 mL·h) −1 , we checked whether other combinations of values of the operational variables falling within the ranges imposed by their errors would lead to maximum values over the range 0.26 ≤ (r A ) global est ≤ 0.27 g acetic acid·(100 mL·h) −1 . A systematic analysis provided the ranges shown in Table 3. With a 0.2 interval for each variable, such ranges led to a total of 162 combinations of which only 109 (see Table S1 in file "S1.docx" in Supplementary Materials) fell within the previous range of (r A ) global est . However, using all 162 combinations ensured that this variable would be very close to its maximum value: 0.25-0.27 g acetic acid·(100 mL·h) −1 . Therefore, any of the 162 combinations would, in theory, lead to the highest possible mean rate of acetic acid formation in the two-bioreactor system-and hence to values of each operational variable falling in the ranges of Table 3. Table 3. Upper and lower end of the ranges over which the operational variables maximized (r A ) global est .  Table 4 (see file "Optimal_rA.m" in Supplementary Materials). The problem was like that stated in Equations (6) and (7) except that the aim was to minimize (r A ) global est rather than maximize it. As in the previous case, the error made in modelling (r A ) global est was used to examine the combinations of values of the operational variables, with provision for their errors (see Table 5), leading to minimum values over the range 0.11 ≤ (r A ) global est ≤ 0.12 g acetic acid·(100 mL·h) −1 . The ranges of Table 5 for the operational variables, with 0.2 unit intervals, led to a total of 550 combinations only 29 of which provided an (r A ) global est value falling in the previous range of minimum values (see Table S2 in file "S1.docx" in Supplementary Materials). However, any of the 550 combinations ensured that (r A ) global est would fall in the range 0.116-0.135 g acetic acid·(100 mL·h) −1 and thus be very close to the minimum of Table 4. Therefore, any combination would provide the lowest mean rate of acetic acid formation in the two-bioreactor system and lead to values of the operational variables falling in the ranges of Table 5. Table 5. Upper and lower end of the ranges over which the operational variables minimized (r A ) global est .  Table 6 for P m est , E u2 est , V u2 est , t cycle est , V m est , EtOH m1 est , EtOH m2 est , HAc m1 est and HAc m2 est under the operating conditions maximizing or minimizing (r A ) global est . The ranges of some variables reflect differences due to the influence of T 2 -like (r A ) global est , other variables were independent of the temperature in the second reactor, however.
H Ac m1 est = 9.188 + 0.0932·E 2 Based on the previous results, the highest and lowest mean rate of acetic acid production were obtained by using extreme values of some operational variables. Thus, the minimum (r A ) global est value (Table 4) was obtained with the lowest temperature in the first bioreactor (T 1 = 28 • C), and the smallest values in the ranges of E l1 and E u1 [viz., 4% and 2% (v/v), respectively]. These results suggest that the smaller the amount of substrate available to AAB and the higher the acidity to which the bacteria are exposed during a fermentation cycle (see EtOH m1 est and HAc m1 est values in Table 6) are, the less suitable will be the medium for their metabolic action -and hence for acetification. Also, because the volume of medium unloaded from the first bioreactor was quite substantial (V u1 = 5.75 L), the amount of inoculum remaining in it for the next cycle was considerably reduced, and so was (r A ) global est as a consequence. This result is consistent with those of modelling studies on a single bioreactor using both first-principles and black-box models [8][9][10]22,23,31,32].
On the other hand, if the aim is to maximize (r A ) global est , the operating conditions should be essentially the opposite of those described in the previous paragraph. As can be seen from Table 6, such conditions resulted in increased mean substrate availability and decreased mean acidity in both bioreactors.
One other interesting conclusion is that P m est is rather different under the conditions maximizing and minimizing (r A ) global est (see Table 6). As expected, the greater (r A ) global est is, the higher will be acetic acid production. However, as shown in the following section, the maximum mean rate did not result in the highest possible P m est level.

Optimization of P m est
The variable (r A ) global est possesses a high industrial interest. In practice, however, the aim of a vinegar producing plant is to obtain as much acetic acid per hour of operation (i.e., to maximize the overall production of acid, P m ). This variable can be estimated with an error of 0.7 g acetic acid·h −1 from Equation (8). In this Equation, only three of the operational variables (E l1 , E u1 and V u1 ) have a direct influence on P m est , but as explained before in the case of Equation (5), not all the terms in Equation (8) have the same influence on P m est (see Table 7); the higher the value of the statistic F, the more significant the term is for P m est ; it can be seen from Table 7 that the most influential term is that of E l1 . Below is described the procedure followed to identify the operating conditions maximizing (the industrial target) or minimizing (an undesirable outcome) P m est .
The procedure of Section 2.3 to develop a MATLAB script to identify the optimum (maximum or minimum) P m est value was used (see file "Optimal_Prod.m" in Supplementary Materials). The conditions maximizing P m est can be identified by solving Equations (17) and (18), the latter containing the constraints imposed by the values of the experimental variables under the experimental conditions used.
where T 2 is the temperature in the second bioreactor. Table 8 shows the value of each operational variable maximizing P m est and hence the solution to the optimization problem.  In a similar way to the previous discussion on (r A ) global est , plotting the response surfaces of P m est varying just a pair of the operational variables at once, it is possible to show that significant changes can only be found when the variable E l1 is included. By way of example, response surfaces varying E l1 and E u1 on one side and E u1 and E l2 on the other are shown in Figure 4a,b, respectively. Values from Table 8 are used for the remaining variables.   Table 8 for the remaining operational variables.
Similarly to ( ) , the errors in and the operational variables were used to identify the combinations of variables leading to values over the range 35.9 ≤ ≤ 36.6 g acetic acid·h −1 (see Table 9). The ranges of Table 9, at 0.2 unit intervals, allowed a total of 4320 combinations to be identified of which 50 provided a value within the previous range (see Table S3 in file "S1.docx" in Supplementary Materials). However, any of the 4320 combinations led to a very high value (32.8-36.6 g acetic acid·h −1 , which are close to that range). Therefore, any of the combinations of operational variables at values within the ranges of Table 9 would lead to a near-maximum value. Minimizing provided the results of Table 10 (see file "Optimal_Prod.m" in Supplementary Materials).  The combinations of values of the operational variables leading to minimum values over the range 14.7 ≤ ≤ 15.4 g acetic acid·h −1 were identified similarly to those maximizing this dependent variable. Using the values in the ranges of Table 11 (0.2 unit intervals) provided 360 combinations of which 36 leading to values within the previous range (see Table S4 in file "S1.docx" in Supplementary Materials). However, all combinations led to very low values: 14.6-17.1 g acetic acid·h −1 , which are very close to the minimum value.  Table 8 for the remaining operational variables.
Similarly to (r A ) global est , the errors in P m est and the operational variables were used to identify the combinations of variables leading to values over the range 35.9 ≤ P m est ≤ 36.6 g acetic acid·h −1 (see Table 9). The ranges of Table 9, at 0.2 unit intervals, allowed a total of 4320 combinations to be identified of which 50 provided a P m est value within the previous range (see Table S3 in file "S1.docx" in Supplementary Materials). However, any of the 4320 combinations led to a very high P m est value (32.8-36.6 g acetic acid·h −1 , which are close to that range). Therefore, any of the combinations of operational variables at values within the ranges of Table 9 would lead to a near-maximum P m est value. Table 9. Upper and lower end of the ranges over which the operational variables maximized P m est . Minimizing P m est provided the results of Table 10 (see file "Optimal_Prod.m" in Supplementary Materials). Table 10. Values of the operational variables needed to minimize P m est . The combinations of values of the operational variables leading to minimum P m est values over the range 14.7 ≤ P m est ≤ 15.4 g acetic acid·h −1 were identified similarly to those maximizing this dependent variable. Using the values in the ranges of Table 11 (0.2 unit intervals) provided 360 combinations of which 36 leading to P m est values within the previous range (see Table S4 in file "S1.docx" in Supplementary Materials). However, all combinations led to very low P m est values: 14.6-17.1 g acetic acid·h −1 , which are very close to the minimum value. Table 11. Upper and lower end of the ranges over which the operational variables minimized P m est .  Table 12 compares the values of the estimated state variables for the maximum and minimum P m est values. As can be seen, the conditions maximizing P m est were different from those maximizing (r A ) global est . In fact, P m est and (r A ) global est are related by the mean volume of medium, V m est , in Equation (19). Since the bioreactors worked in a semi-continuous mode, the mean volume was dependent on the particular operating conditions and those resulting in the maximum possible values of P m est and (r A ) global est need not coincide (see V m est values in Tables 6 and 12).
As stated above for (r A ) global est , an increased mean availability of ethanol and a decreased mean acidity boosted acetic acid production (P m ). However, ethanol and acetic acid concentrations above certain levels have an inhibitory effect whereas too low levels can even boost bacterial activity. These effects have been extensively examined in modelling studies [8][9][10]22,23,31,32]. For example, ethanol concentrations above 6% (v/v) and acetic acid concentrations above 6% (w/v) have been found to have substantial inhibitory effects [8][9][10]. This knowledge, and the mean acidity and ethanol values of Tables 6 and 12, allow the (r A ) global est and P m est results to be more easily understood.
One other major inference is that, based on Table 12, maximizing P m est leaves a substantial amount of unavailable ethanol in the second bioreactor [E u2 = 4.3% (v/v)], which is, in general, industrially unacceptable. It is, therefore, interesting to identify the specific conditions that will maximize P m est while ensuring that the substrate is depleted to an acceptable degree. The procedure followed for this purpose is described in the following section.

Optimizing P m est While Ensuring Enough Substrate Depletion
For economy and legal reasons, the amount of ethanol present in industrially produced vinegar must not exceed prescribed levels depending on each specific type of vinegar. It is therefore important to identify the operating conditions optimizing acetic acid production (P m est in Equation (8)) under a constraint on the ethanol concentration at the time the second bioreactor is unloaded (E u2 est in Equation (9)). The optimization problem is stated in Equations (20) and (21). max P m est s.t. , where constant is the ethanol concentration at unloading time in the second reactor to be imposed. Table 13 shows the values of the operational variables providing the solution to the previous problem with E u2 est values from 0.2% to 1.5% (v/v). Table 13. Values of the operational variables needed to maximize P m est with an E u2 est value of 0.2%, 0.5%, 1.0% or 1.5% (v/v). The solution was obtained by using the MATLAB script of Section 3.2 (see file "Opti-mal_Prod.m" in Supplementary Materials) as expanded with the constraint on E u2 est . The previous range of E u2 est was used to consider virtually complete and less marked depletion of the substrate because some vinegars (particularly those under a designation of origin or aged in casks) are allowed to contain higher levels of residual ethanol.
As expected, the greater E u2 est was (i.e., the less markedly the substrate was depleted in the second bioreactor), the more favored was AAB activity and the higher was acetic acid productivity as a result.

Comparison of the Performance of Serial and Parallel Bioreactors
This was one other point of especial interest. For this purpose, we compared the results obtained with the two bioreactors working in series here and those previously reported for others operating in parallel [8][9][10]22,23,31,32]. Because the previous studies were conducted at 31 • C, we estimated P m and E u2 at that temperature. As can be seen from Table 14, there were no substantial differences in performance between the serial and parallel configurations when P m was maximized under no constraint on E u2 -not even at the optimum T 1 and T 2 values found here. However, maximizing P m with E u2 = 0.5% (v/v) resulted in the serial system being more productive than the parallel one. Specifically, at 31 • C the former system reached a P m 12.2% higher than did the latter. It was impossible to directly compare productivity with the two systems at 32 • C and E u2 = 0.5% (v/v) because no productivity data for the parallel system at that temperature were available. In any case, the results of Table 14 suggest that they cannot be too different from those obtained at 31 • C. Finally, it should be noted that using two serial bioreactors results in very high acetic acid productivity and virtually complete depletion of ethanol. On the other hand, as shown in a previous work [16], depleting the second bioreactor in a system operating in parallel would slow down the process as a result of the time needed for acetic acid bacteria to regain full activity after a period of substrate scarcity and a high acidity in the medium-all of which considerably detract from productivity [16].
One other aim fulfilled by modelling the system and using results of previous research work was to compare the performance of the system with the two reactors operating in series and in parallel. P m can be maximized with no restrictions on E u2 by using either configuration. However, if the substrate is to be depleted [e.g., E u2 = 0.5% (v/v)], then the serial configuration will be more productive.
Supplementary Materials: The following are available online at http://www.mdpi.com/xxx/s1. "Optimal_rA.m" and "Optimal_Prod.m": MATLAB scripts for optimizing (r A ) global est and P m est , respectively. File "S1.pdf": Combinations of operational variables to analyze the maximum and minimum of (r A ) global est and P m est .