Mathematical Modeling and Stability Analysis of a Two-Phase Biosystem

: We propose a new mathematical model describing a biotechnological process of simultaneous production of hydrogen and methane by anaerobic digestion. The process is carried out in two connected continuously stirred bioreactors. The proposed model is developed by adapting and reducing the well known Anaerobic Digester Model No 1 (ADM1). Mathematical analysis of the model is carried out, involving existence and uniqueness of positive and uniformly bounded solutions, computation of equilibrium points, investigation of their local stability with respect to practically important input parameters. Existence of maxima of the input–output static characteristics with respect to hydrogen and methane is established. Numerical simulations using a specially elaborated web-based software environment are presented to demonstrate the dynamic behavior of the model solutions.


Introduction
Anaerobic digestion (AD) is a multi-step biotechnological process widely applied for the treatment of wastes and wastewater, where usually hydrogen (H 2 ) is not accumulated as an intermediary product [1][2][3]. Recently, the interest in H 2 production by modifying operating conditions of the AD has increased. This process has been denominated as dark fermentative H 2 production in contrast with light derived H 2 producing processes [4][5][6][7][8]. In traditional AD, H 2 is not detected as it is consumed immediately e.g., by hydrogenotrophic methanogens to produce methane CH 4 and carbon dioxide CO 2 [9]. Conversely, H 2 can be produced as main component in biogas by engineering the process conditions in a way that favors Clostridium and Enterobacter metabolism against the predominance of archaeal microflora [10,11]. However, the main limitation of dark fermentative H 2 production is the rather low energy recovery, cf. [7,8,12]. In order to completely utilize the organic acids produced during dark fermentation and improve the overall energy conversion efficiency, a two-stage AD (TSAD) concept consisting of a hydrogenic process followed by a methanogenic process has been suggested (see e.g., [12][13][14]).
A lot of models describing separately the fermentative hydrogen production (see, e.g., [15,16]) and the AD for methane production (cf. e.g., [17][18][19][20]) are known. However, only a few models of TSAD processes are available; see, for example, [21][22][23][24] and the references therein. In [25], both bio-hydrogen and bio-methane productions are optimized in two bioreactors regarding two operating parameters, organic loading rate, and dilution rate; however, the first bioreactor is running in semi-continuous and the second-in batch operation mode because the working volumes of both bioreactors are equal. In [26], steady-state analysis for the TSAD process is presented using a very simple model without hydrolysis and without taking into account the most important fact for this scheme-the energy recovery.
In this paper, we propose a mathematical model, describing the process of simultaneous H 2 and CH 4 production by AD of organic wastes in a cascade of two connected continuously stirred bioreactors. The biochemical processes in the first bioreactor include disintegration of organic wastes, hydrolysis, and acidogenesis with hydrogen production. The methane production from acetate (methanogenesis) is separated in the second bioreactor. The proposed model is developed by reducing the well known Anaerobic Digester Model No 1 (ADM1) basic structure elaborated by the International Water Association (IWA) [17]. Modeling of TSAD using ADM1 is also presented in [27], where, however, methane is obtained from both bioreactors.
Our mathematical model is presented for the first time in the conference paper [21]. The present article represents an extension of the latter, providing (i) more detailed assumptions used to construct the model, which are based on practical experiments (Section 3); (ii) more precise and deeper theoretical investigations of the model dynamics such as existence and uniqueness of positive and uniformly bounded solutions, more precise computation of practically important equilibrium points of the model and study of their local asymptotic stability (Section 4); (iii) more detailed investigation concerning the optimization of hydrogen and methane production on steady state operation (Section 5), and showing the existence of maxima, which is important for the applications; (iv) numerical simulations in Section 6, based on a newly developed, so called one-page application "(H2,CH4) Bioreactors ver. 1.2" in the SmoWeb platform.
The main goal of this article is to present rigorous mathematical analysis of the proposed model dynamics, to detect and provide some model-based predictions, which could be useful in designing and engineering real laboratory experiments. The theoretical studies are supported and complimented by a variety of computer simulations.

Process Description
In the TSAD system, relatively fast growing acidogens and H 2 producing microorganisms are developed in the first-stage hydrogenic bioreactor (BR 1 with working volume V 1 , see Figure 1) and are involved in the production of volatile fatty acids (VFA-mainly acetate and butyrate, [8]) and H 2 . Some recent studies on the microbial community in the hydrogen producing BR 1 show that dominant are microorganisms of the genus Clostridium [28]. In BR 1 , pH is kept in the range 5.0-5.5, and acetate degrading methanogens are eliminated. On the other hand, the slow growing acetogens and methanogens are developed in the second-stage methanogenic bioreactor (BR 2 with working volume V 2 , and pH in the range 6.5-8.5), in which the produced VFA are further converted to CH 4 and CO 2 by methane-forming bacteria from the genus Methanobacterium, Methanococcusect, see, for details, [9] and the references therein.
It is known that in the two-stage H 2 + CH 4 system the energy yields are up to 43% more in comparison to the traditional one-stage CH 4 production process (cf. [8,25]).
Assume that the volumes V 1 and V 2 of the bioreactors are constant. Let F 1 and F 2 be the inflows in the first and second bioreactor, respectively, and let F 1 = F 2 = F be valid. It is well known that the dilution rates D 1 and D 2 are defined as Then, (1) Figure 1. Two-phases process of AD with production of hydrogen (H 2 ) and methane (CH 4 ).
Later on in the paper (see Section 5), by analyzing points of maximum biogas production in both bioreactors, biohydrogen in BR 1 and biomethane in BR 2 , we find theoretical values for the hydraulic retention time (HRT) 1/D 1,max and 1/D 2,max , which show that the volume V 2 of the second bioreactor for methane production should be larger than the volume V 1 of the first bioreactor. Therefore, γ < 1 should be valid. We determine the constant γ using the proposed model equations.

Model Description
In this section, we present the mathematical model describing the process of simultaneous H 2 and CH 4 production by AD of organic wastes in a cascade of two continuously stirred bioreactors. The model is derived on the basis of the ADM1 basic structure as well as on our experience with the two-phase AD process with hydrogen and methane production, see [29]. The ADM1 is a complex model, describing all known (till its creation) microbiological and biochemical reactions occurring in AD of organic wastes. ADM1 is suitable mainly for process knowledge and simulation, but not appropriate for theoretical studies and model-based control design due to the plenty of input parameters which are difficult to obtain. The reduction of ADM1 is done here in order to simplify the latter and to provide an opportunity for studying the asymptotic stability of the dynamics with respect to two essential indicators: (i) the separation of the model equations into two groups, corresponding to the reactions taking place in each one of the two bioreactors, and (ii) the type of the organic wastes used.
For simplification of the model, the following assumptions have been accepted: • In BR 1 only the main VFA (acetate) is formed, whereas propionate and butyrate are omitted. In this way, the slow acetogenic phase has been omitted and consequently in BR 2 only methanogenesis occurs.

•
In both bioreactors, the existence of inhibition phenomena is not taken into account.

•
Balance equations of the hydrogen and methane in the liquid phases have been neglected because they are practically not dissolved in liquids.

•
Hydrogenotrophic bacteria, producing methane from hydrogen, are eliminated previously and do not exist in this process.

•
Equations describing balances of inorganic components and some biochemical equations have been neglected in view of simplifying the model. • pH is considered in the range between 5.0 to 5.5 in the first bioreactor BR 1 , and in the interval 6.5-8.5 in the second bioreactor BR 2 . Within these ranges, it is assumed that the negative influence of pH could be disregarded.
These ranges of pH were obtained during our own laboratory experiments with TSAD of waste from maize treatment (not published yet).
The proposed models of the two bioreactors are of balancing type for continuously stirred tank reactors (CSTR) like the original ADM1 model. For each bioreactor, the ADM1 equations are adapted to describe the dynamics of the most important variables for the AD of organic waste.
Following the above assumptions, we present now the mathematical model for hydrogen and methane production in the two bioreactors. The definitions of the phase variables and parameters in the model equations are given in Tables 1 and 2 below. The dynamics in the first bioreactor are described by the following set of 10 nonlinear ordinary differential equations: with gaseous output where Q h2 is the hydrogen flow rate [L/h]. In (2)-(11), the biomass specific growth rates µ suaa,su , µ suaa,aa and µ f a are of the form µ suaa,su = µ suaa,su (S su , S aa ) = k m,suaa S su k s,suaa + S su · S su S su + S aa , µ suaa,aa = µ suaa,aa (S su , S aa ) = k m,suaa S aa k s,suaa + S aa · S aa S su + S aa , We extend the model (2)-(11) by the following two equations describing the process in the second bioreactor (cf. [30]): d dt X ac (t) = (−D 2 + µ ac,ch4 )X ac (13) with gaseous output Q ch4 = Y ch4,ac µ ac,ch4 X ac where Q ch4 is the methane flow rate [L/h], and µ ac,ch4 = µ ac,ch4 (S ac,ch4 ) = k m,ac S ac,ch4 k s,ac + S ac,ch4 .
All model coefficients in the Equations (2)-(13) are assumed to be positive. The dilution rates D 1 and D 2 are the decision or control variables. The numerical values in the rightmost column of Table 2 are taken from [31]. In the model variables and parameters, the subscripts h2 and ch4 indicate hydrogen (H 2 ) and methane (CH 4 ), respectively.

Investigation of the Model Solutions
This section presents a mathematical analysis of the dynamics (2)- (13). First, we establish existence, positivity, uniform boundedness, and uniqueness of the model solutions. The results are presented in Theorem 1. The interested reader can find a detailed proof of Theorem 1 in Appendix A.1 of the Appendix A.
Note that the assumptions 0 < Y suaa < 1 and 0 < Y f a < 1 in Theorem 1 are not restrictive. Since Y suaa and Y f a are yield coefficients, their values are always less than 1 in practice; see last column in Table 2.
In the next two subsections, we find the simplest solutions of the model, namely the equilibrium points, and investigate their local asymptotic stability.

Equilibrium Points of the Model
The equilibrium points of the model are solutions of the algebraic equations obtained from (2)-(13) by setting the right-hand sides equal to zero: We are looking for nonnegative solutions of (14)-(25) since only these equilibrium points correspond to the physical and biological meaning of the model variables. The equilibrium points are computed as functions of the control variables, the dilution rates D 1 and D 2 .
Using Equation (18), we find the equilibrium componentX c with respect to the phase variable X c , which is obviously positive:X Furthermore, Equations (19)- (21) imply the equilibrium componentsX ch ,X pr andX li respectively: Equations (16) and (23) are used to determine the steady state componentsS f a andX f a with respect to the phase variables S f a and X f a . Obviously, Equation (23) suggests that the trivial solution X f a ≡ 0 always exists, but is not of practical interest. Assuming that X f a ≡ 0 and using the expression of µ f a , we obtain from (23) the following steady state componentS f ā It is clear thatS f a exists andS f a > 0 if and only if Equation (16) then delivers the equilibrium component for X f a , Obviously,X f a from (29) satisfiesX f a > 0 if and only if Remark 1. The quantity S in f a + 1 D 1 f f a,li k hyd,liXli on the right-hand side of (30) can be considered as a worst-case upper bound of the total concentration of fatty acids (LCFA) S f a , which accumulation in BR 1 might occur as a result of some imbalances in the degradation phase, leading to wash-out of the LCFA degraders X f a [32].

The inequality (30) is equivalent with
The left hand-side cubic polynomial in (31) always possesses one positive real root D (1) 1 , thenX f a becomes equal to zero; this means that D 1 = D (1) 1 is a (transcritical) bifurcation value of the parameter D 1 , leading to wash-out of the biomass X f a for D 1 > D (1) pay special attention to the caseX f a = 0 motivating our considerations later on in the paper (see Section 5). IfX f a = 0, then we obtain the following steady state component for which always exists and is positive if 1 . Here, and in what follows, the superscript 0 indicatesX f a = 0 in the steady state components. Furthermore, Equations (14), (15) and (22) are used to compute the steady state components S su ,S aa , andX suaa in case positive solutions do exist. It is clear from Equation (22) that the trivial solution X suaa ≡ 0 always exists, but we exclude this case for biological evidence. Let X suaa ≡ 0. Then, Equation (22) reduces to Using Equation (15), we find The equilibrium component with respect to X suaa should be positive and this will be fulfilled if A similar remark to Remark 1 can be made to the above inequality (35) concerning the concentrations of the amino acids S aa and of the degraders X suaa .
Substituting X suaa from (34) into Equation (14) and using Equation (33), we obtain the following nonlinear algebraic system with respect to S aa and S su : It is straightforward to see that system (36) is equivalent to the following one: where All coefficients A i , i = 1, 2, 3, and B j , j = 1, 2, 3, 4, depend on D 1 (the dilution rate in the first bioreactor BR 1 ). Assume that there exists a value D Then, the equilibrium componentX suaa will exist and will be positive for D 1 ∈ 0, min{D Note that an equalityS aa = S in aa + 1 D 1 k hyd,prXpr for some value of D 1 will lead toX suaa = 0 and may cause wash-out of the biomass X suaa , which is practically undesirable, and we shall exclude this case.
Equation (17) is used to find the steady state component The latter is positive if Y suaa < 1 and Y f a < 1 are fulfilled. IfX f a = 0, then the equilibrium component with respect to S ac becomes S 0 ac exists and is positive if Y suaa < 1 and D 1 ∈ D (1) 1 } hold true. The above calculations deliver the following equilibrium points of the model (2)-(11): 1 }.
Consider now Equations (12) and (13) to compute the equilibrium components with respect to X ac and S ac,ch4 . Thereby, we assume thatS ac andS 0 ac are known; see (40) and (41), respectively. Equation (25) has the obvious solution X ac ≡ 0, which is not of practical interest. Assuming that X ac ≡ 0, the equality µ ac,ch4 (S ac,ch4 ) = D 2 delivers the following steady state component for S ac,ch4 Obviously, there existsS ac,ch4 > 0 if D 2 < k m,ac holds true. Then, the corresponding equilibrium componentX ac isX The latter is positive ifS If we considerS 0 ac instead ofS ac we obtain which is positive ifS 0 ac >S ac,ch4 .
Therefore, the extended model (2)-(13) can possess the following two equilibrium points, parameterized on D 1 and D 2 : We summarize the above calculations in the following theorem.
Theorem 2 can be specified using the coefficient values in Table 2. Numerical calculations produce the following values: The left plot in Figure 2 visualizes the equilibrium componentS aa , which exists for D 1 ∈ (0, D 1 ), as well as the upper bound S in aa + 1 D 1 k hyd,prXpr ofS aa according to (35). Obviously, (35) is satisfied 1 ) so that D 1 . The right plot in Figure 2 shows the equilibrium component X suaa (D 1 ), which vanishes at D 1 = D (2) 1 . 2 . This implies the following relation: 2 ) = (0, 0.0156).
The equilibrium E(D 1 , D 2 ) corresponds to coexistence of all microorganisms in both bioreactors BR 1 and BR 2 , whereas E 0 (D 1 , D 2 ) is related to wash-out of the LCFA degrader X f a in BR 1 for hydrogen production.

Local Asymptotic Stability of the Equilibrium Points
We shall first investigate the local asymptotic stability of the coexistence equilibrium E 1 (D 1 ) for D 1 ∈ 0, D (1) 1 and of the wash-out steady state (withinX f a = 0) E 0 1 , i.e., of the equilibrium points corresponding to the model (2)-(11) of the first bioreactor BR 1 for hydrogen production.
It is well known that an equilibrium point is locally asymptotically stable, if all eigenvalues of the Jacobi matrix evaluated at this equilibrium have negative real parts, cf. e.g., [33]. The eigenvalues of this Jacobi matrix coincide with the roots of the corresponding characteristic polynomial.
Detailed calculations of the characteristic polynomials P(E 1 ; λ) and P(E 0 1 ; λ) with respect to the equilibrium points E 1 = E 1 (D 1 ) and E 0 1 = E 0 1 (D 1 ) are given in Appendix A.2 of the Appendix A. We present them here for reader's convenience: 1 ), where R(E 1 ; λ) is the characteristic polynomial of the sub-matrix ∆ (3) evaluated at the steady state 1 ), and R(E 0 1 ; λ) is the characteristic polynomial of the sub-matrix ∆ (3) evaluated at the steady state components (S su ,S aa ,X suaa ) for D 1 ∈ (D (1) The linear quantities in P(E 1 ; λ) = 0 and P(E 0 1 ; λ) = 0 produce negative real eigenvalues of the Jacobi matrices. In the first polynomial P(E 1 ; λ), the quadratic equation It remains to show that R(E 1 ; λ) = 0 and R(E 0 1 ; λ) = 0 have roots with negative real parts. These two problems are solved numerically using the coefficient values in Table 2. The real parts of the eigenvalues are presented graphically in Figures 4 and 5, and are obviously negative. Therefore, the steady states 1 ), are locally asymptotically stable.
The positive signs of A, B and A 0 , B 0 imply that the above quadratics possess two roots with negative real parts. This proves the local asymptotic stability of the equilibrium points E(D 1 , D 2 ) for D 1 ∈ (0, D 2 ) and E 0 (D 1 , D 2 ) for D 1 ∈ (D (1)

Optimization of Hydrogen and Methane Flow Rates
In this section, we shall find values of the dilution rates D 1 and D 2 , so that the hydrogen and methane flow rates Q h2 and Q ch4 evaluated at the model equilibrium points achieve maximal values.
We compute Q h2 on the two sets of steady states E 1 (D 1 ), D 1 ∈ (0, D (1) Recall that the superscript 0 in the equilibrium E 0 1 (D 1 ) indicatesX f a = 0. The functions Q h2 (D 1 ) and Q 0 h2 (D 1 ) depend on the decision (manipulated) input D 1 and are called input-output static characteristics with respect to the hydrogen production.
The important (from practical point of view) question is whether the functions Q h2 (D 1 ) and Q 0 h2 (D 1 ) are unimodal with respect to D 1 in the admissible intervals for D 1 . For example, the function Q h2 (D 1 ) is called unimodal if there exists a unique (admissible) point D 1,max , such that Q h2 (D 1 ) possesses maximum Q h2,max = Q h2 (D 1,max ), Q h2 (D 1 ) is strongly increasing if D 1 < D 1,max and Q h2 (D 1 ) is strongly decreasing if D 1 > D 1,max . Figure 6 (left plot) presents the graphs of the input-output static characteristics Q h2 (D 1 ) for 1 ) coincides with the maximum of Q 0 h2 (D 1 ). From a practical point of view, this means that the maximal flow rate of H 2 could be achieved at a operation mode, related to wash-out of the LCFA degraders X f a in BR 1 .
2 ) (right). The vertical dash-dot lines in the left plot pass through D  A similar question of maximizing the biogas production at steady states has been discussed in [34] for the well known four-dimensional AM2 model. There, two maxima of the input-output static characteristics with respect to methane flow rate have also been determined, related to coexistence or wash-out of the acidogenic bacteria in the bioreactor.
Using the methane flow rate Q ch4 , we compute further the input-output static characteristics Q ch4 and Q 0 ch4 on the set of the steady states E(D 1 , D 2 ) and E 0 (D 1 , D 2 ) respectively, namely: The latter are parameterized on the two inputs D 1 and D 2 , sinceX ac andX 0 ac depend on both D 1 and D 2 ; see (43) and (45), respectively.
With D 1 = D 0 1,max ≈ 0.0439, we compute the equilibrium componentsS ac,ch4 andX 0 ac according to (42) and (45), as well as Q 0 ch4 , all of them as functions of D 2 . Figure 6 (right plot) visualizes the graph of Q 0 2 ); the latter is a unimodal function, taking its maximum at D 0 The latter results are theoretical and predicted by the model studies. Experimental data that confirm these results by considering H 2 and CH 4 volumetric production under similar conditions for which our model was developed (i.e., two bioreactors in a cascade with a similar ratio of working volumes and operating in continuous mode) can be found in [35].
Using the presentation (1), we define and compute the constant γ and in this way the relationship between the volumes V 1 and V 2 of the two bioreactors: The last relation implies that the volume of the second bioreactor BR 2 for methane production should be at least 4.4 times greater than the volume of the first bioreactor BR 1 for hydrogen production.

Dynamic Behavior of the Model Solutions
The dynamic model has been implemented as one-page application "(H2,CH4) Bioreactors ver. 1.   Figure 7 presents the time evolution of the flow rates Q h2 (t) and Q ch4 (t) when D 1 and D 2 take values from Table 3. Figures 8 and 9 display the model solutions by variable D 1 and D 2 given in Table 3. It is seen that, if D 1 ∈ (0, D     Table 3.  Table 4. The model solutions tend to the corresponding components of the coexistence equilibrium E(D 1 , D 2 ), which is the locally stable steady state. These numerical simulations demonstrate the fact that uncertainties in the input X in c do not affect the stability of the model dynamics.    Table 4.

Discussion
In this paper, we propose a mathematical model describing the process of simultaneous production of hydrogen and methane by anaerobic digestion of organic wastes in a cascade of two connected continuously stirred bioreactors with different volumes. The proposed model is developed by adapting and reducing the universal Anaerobic Digester Model No 1 (ADM1). The equations in our model are separated in two groups, corresponding to the processes in the two bioreactors, for hydrogen and for methane production, respectively. The simplifications made in the two models naturally lead to some limitations compared to the original ADM1. However, as already shown, the latter also does not cover all possible reactions, e.g., possible inclusion in the model of the activity of the so called syntrophicacetate oxidizers [36]. Some new, more complex models have also been developed [37], which however have the same main disadvantage as ADM1-great complexity and thus impossibility to carry out some mathematical analytical studies.
The investigation of the hydrogen and methane flow rates, Q h2 and Q ch4 respectively, on steady state operation, shows existence of maxima for some values of the dilution rates D 1 and D 2 in the two bioreactors BR 1 and BR 2 . The local maximum for D 1 in BR 1 reflects the maximal yield of hydrogen from lipids at very long hydraulic retention time HRT = 1/D 1,max ≈ 1/0.0072 ≈ 139 [hours] ≈ 5.8 [days]. This maximum is much less than the maximum yield of hydrogen derived from monosaccharides and aminoacids (proteins) for the accepted combination of coefficients and initial conditions reflecting the biodegradation of organic wastes. In [35] under similar conditions, results close to our theoretical model-based predictions are obtained.
The model also allows us to find the optimal ratio between the volumes V 1 and V 2 of the two bioreactors subject to the same optimization goal. Many studies in the literature report on TSAD processes; however, most of them are operated manually and the ratio of working volumes of CSTRs is not discussed. In some references, a ratio of 10 is accepted without comments. In [38], a TSAD system consisting of two CSTRs operating at mesophilic conditions are used to investigate the effect of HRT on hydrogen and methane production; the ratio of working volumes of the bioreactors is equal to 8 (without comments). In [39], optimization of TSAD of separately collected municipality bio-waste is carried out in pilot scale CSTRs at thermophilic temperature, using recirculation of the digestate in the second bioreactor to maintain pH in the first bioreactor at optimal value. There, the ratio of working volumes of the bioreactors is equal to 3.8 (without comments). In [40], optimal loading rate is obtained providing maximum H 2 and CH 4 productions in TSAD of cassava wastewaters using specific thermophilic bioreactors and constant recycling ratio of 1:1 with automatic control of pH in the hydrogenic bioreactor. The ratio of working volumes of both bioreactors is equal to 6 (without explanations). The review article [41] is the first one that combines optimization approaches for three possible products from AD-methane, hydrogen and volatile fatty acids (VFA), taking into account different process parameters and types of bioreactors, including acidogenesis and methanogenesis separation in two different bioreactors. However, the ratio of the bioreactors working volumes is not given and discussed.
Using the adapted values of the coefficients in our model, we obtain an optimal ratio of the working volumes of both bioreactors equal to 4.42 using a criterion for maximizing the bioenergy production. This value falls into the range of the above reported values; however, it is theoretically derived.

Conclusions
In this paper, we consider a mathematical model describing a biotechnological process of hydrogen and methane production by anaerobic digestion. The process is carried out in two connected continuously stirred bioreactors, thus the model equations are separated in two groups, describing the processes in each bioreactor. We establish existence and uniqueness of positive and uniformly bounded solutions of the model dynamics. Two equilibrium points of the model are computed and their local asymptotic stability is studied with respect to the practically important input parameters, the dilution rates D 1 and D 2 in the two bioreactors. The first equilibrium point corresponds to coexistence of substrates and biomass, and the second one is related to wash-out of the LCFA degraders in the first bioreactor. We show that there exists a value of the dilution rate D 1 , at which the coexistence equilibrium disappears as a result of bifurcation of the steady state and the wash-out equilibrium remains the only locally asymptotically stable equilibrium. The investigation of the hydrogen and methane flow rates, Q h2 and Q ch4 , on steady state operation, shows existence of maxima at some values of the dilution rates D 1 and D 2 . It is worth noting that the maximal hydrogen flow rate is achieved at the steady state where the LCFA degrader X f a is washed out. This is established here by taking numerical values for the model parameters from the literature [31]. Nevertheless, this fact could help in the design of an AD process to achieve maximal production of either hydrogen and methane. The future work envisages extending the model of the second bioreactor by equations describing an acetogenesis phase, where other VFA (propionate, butyrate and valeriate) will be transformed into acetate, as well as including inhibitory kinetics of the microorganisms growth. The elaborated web-based software will allow us to check different model-based optimization and control strategies that will contribute to designing and engineering a real process. Equation (6) is linear with respect to X c and can be solved explicitly: The solution X c (t) obviously exists for all t ≥ 0, is positive, and is uniformly bounded. In practice, it is reasonable to consider only positive initial conditions X c (0) > 0.
Since X ch (t) > 0 for all t ≥ 0, it follows from Equations (6) and (7) thaṫ Multiplying both sides of the inequality by e D 1 t > 0, we obtain e D 1 tΣ The latter inequality means that lim sup t→∞ Σ 1 (t) ≤ Σ in 1 . Since X c (t) is positive and bounded, and X ch (t) is positive, this implies that X ch (t) is uniformly bounded and exists for all t ∈ [0, ∞).
Denote Σ 2 = X pr + f pr,xc X c , Σ in 2 = f pr,xc X in c . Then, using the fact that X ch (t) ≥ 0 for all t ≥ 0, we obtain from Equations (6) and (9) which means that lim sup t→∞ Σ 2 (t) ≤ Σ in 2 . Since X c (t) is positive and bounded, the latter inequality implies that X pr (t) is also uniformly bounded, and exists for all t ∈ [0, ∞).
Define Σ 3 = X li + f li,xc X c and Σ in 3 = f li,xc X in c . Then, using Equations (6) and (9) and the fact that X li (t) ≥ 0 for all t ≥ 0, we obtaiṅ which implies lim sup t→∞ Σ 3 (t) ≤ Σ in 3 . Since X c (t) is positive and bounded, the latter inequality means that X li (t) is also uniformly bounded, and exists for all t ∈ [0, ∞).
Consider now Then, we obtain from Equations (6), (9), and (11) Since X c (t) and X li (t) are positive and bounded, S f a (t) and X f a (t) are positive, the latter inequality implies that S f a (t) and X f a (t) are uniformly bounded and exists for all t ∈ [0, ∞).
Further consider Σ 5 = S su + S aa + 1 Y suaa X suaa + X ch + X pr + f su,li X li + f ch,xc k dis + f li,xc f su,li + f pr,xc X c , Σ in 5 = S in su + S aa,in + f ch,xc k dis + f li,xc f su,li + f pr,xc X in c .
Then, we obtaiṅ Since X c (t), X ch (t), X li (t) and X pr (t) are positive and bounded, the above presentation implies that S su (t), S aa (t) and X suaa (t) are also uniformly bounded, and exist for all t ∈ [0, ∞).
Consider now Equation (5). Since S ac (t) serves as input into the second bioreactor, see Equation (12), it is reasonable to consider only positive initial values S ac (0) > 0. The positivity of the phase variables included on the right-hand side of Equation (5) imply thatṠ ac (t) ≥ −D 1 S ac (t) and thus S ac (t) ≥ S ac (0)e −D 1 t > 0 for all t ∈ [0, ∞).
Furthermore, the boundedness of the specific growth rate functions imply that there exists a constant Γ 1 > 0 such thatṠ ac (t) ≤ −D 1 S ac (t) + Γ 1 for each t ≥ 0. Then, we have consecutively which means that lim sup t→∞ S ac (t) ≤ Γ 1 D 1 .
Therefore, all solutions of the model (2)-(11) are positive, uniformly bounded, and thus exist for all t ≥ 0.
. . . Using the explicit dependance of g j , j = 1, 2, . . . , 10, on the above variables, we obtain The characteristic polynomial corresponding to the Jacobi matrix J is defined by det(J − λI 10 ), where λ is any complex number and I 10 is the (10 × 10)-identity matrix.