Mathematical study of a two-stage anaerobic model when the hydrolysis is the limiting step

a


Introduction
Anaerobic digestion is the degradation of organic matter in the absence of oxygen.The final product being methane, a renewable energy.This process is more and more used for the treatment of liquid and solid waste.Because of its relative instability due to the possible accumulation of intermediate products, notably the volatile fatty acids (VFA), the modeling of this process has been extensively studied over these last years.Such models are multi-step mass balance models in which the reactional network consists in a number of biological reactions taking place the medium in parallel and / or in series.Their complexity highly depends on the objectives by the modeler.On the one hand, when the objective is to develop models for integrating and formalize the available knowledge typically to better understand bioprocesses, models are generally high order models and not tractable from a mathematical viewpoint, cf. for instance the ADM1 [1].On the other hand, when the aim of the modeling is to develop decision tools or control systems, low order models are better suited, as for instance the AM2 [3].In this second class of models, several two-steps models have been proposed in the literature.Two steps models are commonly used to describe commensalistic microbial systems which take the form of a cascade of two biological reactions where one substrate S 1 is consumed by one microorganism/ecosystem X 1 to produce S 2 which serves as the main limiting substrate for a second microorganism/ecosystem X 2 as schematically represented by the following reaction scheme representing a simplified scheme of the anaerobic digestion The most general two-step model under interest in the actual paper can be written as : where D is the dilution rate, while S in 1 and S in 2 are the input substrate concentrations respectively.Parameters Y i are yield coefficients associated to the bioreactions, k i are the mortality terms while α ∈ [0, 1] is a term allowing to decouple the retention time applying on substrates (supposed to be soluble) and biomass (supposed to be particulate).The kinetics µ 1 and µ 2 are of Contois and Haldane type.By Contois or Haldane-"type", we mean functions that are defined by qualitative properties.For µ 1 , it means a density dependant function which is increasing with S 1 but decreasing with X 1 while µ 2 is non monotonic (cf.modeling hypotheses in the next section).The different analyses of the class of models (1) available in the literature essentially differ on the way the growth rate functions are characterized and whether a specific input for S 2 is considered or not (i.e., the presence of a term S in 2 in the dynamic equation of S 2 ).They differ also in the presence or not i) of a coefficient α allowing to decouple the solid and liquid retention times and ii) of a maintenance term k i (.).For details and informations on the various models considered in the existing literature the reader can refer to [9], or Table (2) and Table (3) in the review paper [13].
A first mathematical study of a pure commensalistic model was proposed by Reilly [7], where µ 1 : and some extensions.He was interested at explaining surprising oscillations observed during an experiment realized in making Saccharomyces carlsbergensis growing on fructose produced by Acetobacter suboxyduns from mannitol.In particular, he established theoretical conditions involving an interval feedback from the yeast to the bacteria.For more information on commensalim, the reader can refer to Stephanopoulos [10].An important contribution on the modeling of anaerobic digestion as a commensalistic system is the model by Bernard et al. [3].The authors considered a Monod function for µ 1 and a Haldane function for µ 2 .Sbarciog and al. [8] studied this model for α = 1 while the most interesting case where 0 < α < 1 and where growth functions were characterized by qualitative properties was studied in [2].
Depending on the nature of treated waste, the limiting step of anaerobic digestion is not the same.If treated waste is liquid, the main limiting step usually considered to be the methanogenesis: in such a case, simple models including only acidogenesis and methanogenesis can be used (cf.for instance Bernard et al. [3]).If the waste contains a high proportion of solid matter, it is the rule rather than the exception to consider that the hydrolysis which is the main limiting step and a model including only hydrolysis and methanogenesis can be used.However, as proposed in Vavilin, [11], when hydrolysis is the limiting step, rates depending only on substrate concentrations, such as Monod model, are not the most appropriate.It is better to describe such complex phenomena by density-dependent models, such as density dependent kinetics, a family in which Contois models falls.Contois model.Using this model, the rate of the hydrolysis step is modeled as : which exhibits the following properties specific to hydrolysis (cf.Mottet et al. [6]) While the analysis of the general model of anaerobic digestion initially purposed by Bernard et al., [3] (representing acidogenis and methanogenesis steps) has been realized be Benyahia et al., [2], from the best of authors knowledge, a two-step model where the kinetic of the first step is modeled by a density-dependent kinetics while the second step exhibits a "Haldane-type" function has never been studied in the literature.It is the aim of the actual paper to study such a generic model.
The paper is organized as follows.In section 2 we present the two-step model with two input substrate concentrations and we give the general hypotheses on the growth functions.In section 3 we give the expressions of the steady states and in section 4, we discuss their stability.In section 5, we illustrate the effect of the second input substrate concentration on the steady states, in designing the operating diagrams, first, with respect to the first input substrate concentration and the dilution rate and second, with respect to the second input substrate concentration and the dilution rate.

Mathematical model
The two-step model reads: where S 1 and S 2 are the substrate concentrations introduced in the chemostat with input concentrations S in 1 and S in 2 .D 1 = αD + k 1 and D 2 = αD + k 2 , where D is the dilution rate, k 1 and k 2 represent maintenance terms and parameter α ∈ [0, 1] represents the fraction of the biomass affected by the dilution rate while Y i are the yield coefficients.X 1 and X 2 are the hydrolytic bacteria and methanogenic bacteria concentrations.The functions are the specific growth rates of the bacteria.
As underlined in the introduction, particular kinetics models as Contois function verifies H1 while Haldane function verifies H2.
A2. f 2 (s 2 ) is positive for s 2 > 0, and satisfies f 2 (0) = 0 and f 2 (+∞) = 0. Moreover f 2 (s 2 ) increases until a concentration s M 2 and then decreases, with Model (4) has a cascade structure which renders its analysis easier.In other terms s 1 and x 1 are not influenced by variables s 2 and x 2 and their dynamics is given by: The behavior of this system is well-know, see [4].A steady state (s * 1 , x * 1 ) must be solution of the system From the second equation we deduce that x * 1 = 0, which corresponds to the washout E 0 = (s in 1 , 0), or s * 1 and x * 1 must satisfy both equations Let γ a function defined by : According to the hypothesis A1, γ(s 1 ) is strictly increasing over the interval ]0, s in  Hence, for x * 1 = 0, the equilibrium E 1 (s * 1 , x * 1 ) exists if and only if D 1 < f 1 (s in 1 , 0).The local stability of the steady state is given by the sign of the real part of eigenvalues of the Jacobian matrix.In the following, we use the abbreviations LES for locally exponentially stable.
Proposition 2.1.Assume that assumptions A1 and A2 hold.Then, the local stability of steady states of ( 5) is given by : When E 0 and E 1 coincide, the equilibrium is attractive (the eigenvalue equal to zero) The results are summarized in the following table : Stable when it exists Proof The local stability of (5) of each steady state depends on the sign of the real parts of the eigenvalues of the corresponding Jacobian matrix.At a given steady state (s 1 , x 1 ), this matrix is given by : where The eigenvalues of J are the roots of its characteristic polynomial det(J − λI).
• For E 0 = (s in 1 , 0), the Jacobian matrix reads Its eigenvalues are λ 1 = −D and For being stable, we need λ 2 < 0. Therefore, E 0 is stable if and only if • For E 1 = (s 1 , x 1 ), where f 1 (s 1 , x 1 ) = D 1 , the Jacobian matrix becomes : As the eigenvalues λ 1 and λ 2 of a square matrix of dimension two are the solutions of the equation The real parts are strictly negative if and only if detA > 0 and T rA < 0.
Therefore by A1 we have

Operating diagram
Apart from the two operating (or control) parameters, which are the input substrate concentration s in 1 and the dilution rate D that can vary, all others parameters (α, k 1 and the parameters of the growth function f 1 (s 1 , x 1 )) have biological meaning and are fixed depending on the organisms and substrate considered.The operating diagram shows how the system behaves when we vary the two control parameters s in . Therefore, the curve separates the operating plan in two regions as defined in Fig. (2).The curve Γ is the border which makes E 0 unstable and at the same time E 1 exists.The following table indicates the stability properties of steady state in each region.
Except for small values of D and s in 1 , notice that the operating diagram of this first part of the 2-step system under study is qualitatively similar to that one of the first part of the AM2 model (when a Monod-like growth function is considered, cf. [3]).

The dynamics of s 2 and x 2
Due to the cascade structure of (4), the dynamics of the state variables s 2 and x 2 are given by with where s 1 (t), x 1 (t) is a solution of (5).A steady state (s * 1 , x * 1 , s * 2 , x * 2 ) of ( 4) corresponds to a steady state (s * 2 , x * 2 ) of ( 10) where either (s where The first case corresponds to (s * 1 , x * 1 ) = E 0 and the second to The system (11) corresponds to a classical chemostat model with Haldane-type kinetics, including a mortality term for x 2 and an input substrate concentration depending on the input flow rate.The behavior of such a system is well-know, see [4].A steady state (s * 2 , x * 2 ) must be a solution of the system From the second equation we deduce that x * 2 = 0, which correspond to the washout In general, by the hypothesis A2, and if we suppose this equation has two solutions such that s 1 2 < s 2 2 .Therefore the system has two positive steady states F 1 = (s 1  2 , x 1 * 2 ) and F 2 = (s 2 2 , x 2 * 2 ), where For i = 1, 2, the steady states F i exist if and only if s in * 2 > s i 2 .Proposition 2.2.Assume that assumptions A1, A2 and (15) hold.Then, the local stability of steady states of ( 11) is given by : 2 ).Proof The Jacobian matrix of the system (11) is written as follows • For F 0 = (s in * 2 , 0), the Jacobian matrix reads The eigenvalues of J F0 are λ 1 = −D and • For F 1 = (s 1 2 , x 1 * 2 ), the Jacobian matrix reads 2 ), the Jacobian matrix reads 2 is negative, hence the equilibrium F 2 is unstable because the eigenvalues are of opposite signs.This completes the proof of the proposition.
The results are summarized in the following table

Operating diagram
The operating diagram shows how the system behaves when we vary the two control parameters s in * . Therefore, the horizontal line together with the curve separate the operating diagram plane in three regions as defined in Fig.The following table indicates the stability properties of steady states in each region.

Steady state of (4)
The aim of this section is to study the dependance of the steady state of (4) with respect to the operating parameters D, s 2 ), and we denote 2 ), and we denote 2 ), and we denote 2 ), and we denote 2 ).We summarize the results obtained in the previous sections in the following proposition.
Proposition 2.3.The system (4) has at most six steady states : 2 ), exists if and only if s in 2 > s 1 2 .
2 ), exists if and only if s in 2 > s 2 2 .

Stability of steady states of (4)
Now we will study the stability of the steady states given in the Proposition (2.3).For this we consider the Jacobian matrix where J 11 is defined in (8) and J 22 is defined in (17).This matrix has a blocktriangular structure.Hence, the eigenvalues of J are the eigenvalues of J 11 and the eigenvalues of J 22 .The existence of steady states depend only on the relative positions of the two numbers s in 1 and s * 1 defined by ( 7) and of the four numbers s 1 2 , s 2 2 defined by ( 14), s in 2 and s in * 2 solutions of (12).Equilibria stability is given in the following table. ).
T r(J 11 ) < 0 and det(J T r(J 11 ) < 0 and det(J 11 ) > 0 by A1 T r(J 11 ) < 0 and det(J 11 ) > 0 These results are summarized in the following tables, where S and U read for LES and unstable respectively and no letter means that the steady state does not exist.
Remarque 2.4.Here we have excluded the limit values in the stability condition (Ex: ..) These cases are related to the eigenvalues of the Jacobian matrix with a real part equal to 0 and the corresponding states are named nonhyperbolic stationary states.Otherwise, they are named hyperbolic steady states.
The different possible cases of non-hyperbolic equilibria are summarized in the following theorem Theorem 2.5.If s in 1 < s * 1 , then we have three subcases Let us give the details of the proof in the case 2.9.The other cases can be studied similarly.Assume that s in 2 = s 1 2 < s in * 2 < s2 2 , then x 1 2 = 0, x * 1 > 0 and x 1 * 2 > 0 by (16).Therefore (see Proposition 2.3), the system has three equilibria 2 and E 1 2 .Using the linearization, we obtain that E 0 2 and E 1 2 are hyperbolic.
Remarque 2.6.In each case among those cited in the table (2), ( 3) and the previous Theorem 2.5 we associate a corresponding figure which represents the relative position of the following parameters s 1 2 , s 2 2 , s in 2 and s in *

Simulations
To illustrate our results, we plot what is called the 'operating diagrams' in a number of situations for the system (3) under hypothesis H1 and H2.Recall that operating diagram summarizes the existence and the nature of the equilibria of a dynamical system as a function of its input variables.Here, these control inputs are D, S in 1 and S in 2 .More particularly, we plot either the operating diagrams in the plan {S in 1 , D} for a fixed value of S in 2 or in the plan {S in 2 , D} for a fixed value of S in 1 .All simulations are performed with the following growth functions: The choice of the values of the model parameters is difficult.Here, our objective is not to match a set of data but rather to highlight the most interesting qualitative properties of the model under interest.To do so, most parameters are taken from [3], while others were changed more significantly as I. Indeed, as underlined in [5], inhibition of the second reaction is not visible if original parameters proposed in [3] are used considering reasonable ranges of variations for S 1 and S 2 .With respect to this later, the Haldane parameter I was thus significantly decreased to willingly increase the effect of the inhibition of the growth of X 2 by S 2 .Finally, parameter values used are summarized in the following table To compute the different regions of the operating diagrams, we use the numerical method reported in [5] which the algorithm is recalled hereafter.

Algorithm for the determination of the operating diagrams
The algorithm is as follows: for each value of input variables chosen on a grid, the equilibria are computed.The eigenvalues of the Jacobian matrix are then calculated for each equilibrium.Finally, according to the conditions of existence and the sign of the real parts of the eigenvalues, a 'flag' is assigned to each of the 6 equilibria: 'S' for stable, 'U' for unstable or nothing if it does not exist.This procedure stops when all the values of the grid {S in , D} have been scanned.As a result, we obtain a number of 'signatures' composed of sequences of 'S', 'U' or 'nothing' coding for the existence and stability of the equilibria that we group into regions as summarized in the tables at the end of sections 2.1 and 2.2, respectively.This algorithm may be formalized as follows: let N 1 , N 2 be two integers in N * and h 1 = D N1 and h 2 = Sin N2 the two iteration steps:

Operating diagrams
In this part we illustrate our results by plotting the operating diagrams and discuss them. 2 and E 1 2 , and A 7 (in dark blue) is the stability region of steady-state E 1 2 , the difference between the regions A 5 and A 7 being that the equilibrium E 1 1 does not exist in the region A 5 but exists and is unstable in A 7 , see tables (2) and (3).This last diagram allows us to see the appearance/disappearance of steady states as a function of the input variable D (recall that S in 1 and S in 2 are fixed).As long as D is small enough (i.e.such that αD + k 2 < d 1 ), the quantity of substrate entering the second step of the reaction is very important: the system is in the region A 7 where the positive equilibrium is the only stable equilibrium.As D increases the size of the attraction basin of this equilibrium decreases until D reaches a critical value (corresponding to the point P 1 in Fig. (7), the frontier between regions A 7 and A 6 ).This critical value corresponds to that one for which the term S in * 2 becomes exactly the largest solution of the equation µ 2 (S 2 ) = D 2 (equivalent to equation (14) for the system (4)): the system enters then in the region A 6 .From a biological point of view, the interpretation is as follows: as D increases, X * 1 decreases and thus S in * 2 decreases as can be seen from equation (12).When D 2 = d 2 , the quantity of available resource necessary to X 2 to grow may become limiting for some initial conditions, leading the system to enter a bistability zone.With the values of the parameters chosen, further increasing D leads definitely X 2 to the washout: the system enters into A 5 in crossing the point P 2 of Fig. (7).Finally, if D is such that D 1 = d 2 (the critical value corresponding to the maximum growth rate of X 1 ) X 1 goes also to extinction and the system enters into A 1 .
Example 2: Let S in 1 = 14 g/L.This case is even more interesting since, when D inscreases, the system goes back to A 5 once before leaving it definitely in browsing the following regions: A 7 A 6 A 7 A 5 A 1 .While D is small enough (i.e.such that αD+k 2 < d 3 , cf.Fig. ( 9)), the reasonning remains the same than before: the only difference is that the value of D leading the system to enter into A 6 through P 3 is a little bit higher than in the previous case (d 3 > d 1 ).It is due to the fact that the second step of the reaction receives less input from the first step when compared to the case where S in 1 = 18g/L, thus enlarging the attraction basin of the stable positive equilibrium.Then, when D is further increased, it may happen an interesting phenomenon: the system enters back into A 7 through point P 4 instead of entering A 5 as it was the case before.In fact, this strongly depends on model parameters and in particular on the relative rate at which S in * 2 and the largest solution of the equation (12) vary as functions of D see Fig. (10).In other words, it depends on how the input concentration of the second step S in * 2 -which includes the part of S 1 transformed into S 2 during the first reaction -is affected by D. On the one hand, if we are in a 'flat' zone of the growth rate µ 2 assuming we consider concentrations at the right of the maximum of µ 2 , a small variation of D (and thus of D 2 ) will change very much the largest solution of the equation (12) while S in * 2 may almost remain constant.On the other hand, if D is such that D 2 crosses µ 2 in a sharper zone of the Haldane function (typically around the inflexion point, still considering S 2 evolves at concentrations such as we are on the right of the maximum of µ 2 ), a small change on D (and thus on D 2 ) will affect much more the solution of the equation ( 12) than before.In any of these situations, the relative positions of the largest solution of the equation ( 12) and of S in * 2 determine whether the system will evolve in the region A 7 or A 6 and we can observe, as D increases to a value such αD + k 2 = d 4 a return of the system from A 7 into A 6 through P 3 and then back into A 7 through P 4 .It is well illustrated in the bifurcation diagram plotted in Fig. (10).

1
and D. The operating diagram is shown in Fig. (2).The condition

2 and
D. The operating diagram is shown in Fig. (3).The conditions s in * 2

Fig. ( 4 )
Fig. (4) represents the operating diagram of model (1) in the plan {S in 1 , D} for a fixed value of S in 2 .The regions are defined as follows.A 1 (in green) is the stability region of the washout E 0 1 , A 5 (in blue) is the stability region of steady-state E 1 2 , A 6 (in red) is the bistability region of the steady-states E 02 and E 1 2 , and A 7 (in dark blue) is the stability region of steady-state E 1 2 , the difference between the regions A 5 and A 7 being that the equilibrium E1  1 does not exist in the region A 5 but exists and is unstable in A 7 , see tables (2) and (3).

Fig. ( 5 )Figure 6 :
Fig. (5) represents the operating diagram of model (1) in the plan {S in 2 , D} for a fixed value of S in 1 .The regions are defined as follows.A 2 (in yellow) is the stability region of steady-state E 1 1 , A 3 (in orange) is the bistability region of the washout E 0 1 and the steady-state E 1 1 , A 8 and A 9 (in pink) and (in dark pink), respectively,

Figure 8 :
Figure 8: The Bifurcation diagram for the input control D for S in 1 = 18g/L.

Table 1 :
The stability conditions for the system (4)

Table 2 :
The three cases when s in 1

Table 4 :
Parameters values of the simulation Algorithm 1 Operating diagram for i varying from 1 to N 1 do; for j varying from 1 to N 2 ; determine 6 equilibria of the model E 1 ...E 6 for k varying from 1 to 6 do calculate the jacobian matrix at E k (J E k ) calculate the eigenvalues of (J E k ) if all (the conditions of existence of E k are fulfilled and all real parts of the eigenvalues of (J E k ) are non-positive) then E k is stable else if all (conditions of existence of E k are fulfilled at least one real part of eigenvalue of (J E k ) is positive) then E k is unstable else E k does not exist