{'fraction_non_alphanumeric': 0.04501757855416392, 'fraction_numerical': 0.03425071412876291, 'mean_word_length': 3.7574807265124788, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 38, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Operating Diagram for a Two-Step Anaerobic Digestion Model

The Anaerobic Digestion Model No. 1 (ADM1) is a complex model which is widely accepted as a common platform for anaerobic process modeling and simulation. However, it has a large number of parameters and states that hinder its analytical study. Here, we consider the two-step reduced model of anaerobic digestion (AM2) which is a four-dimensional system of ordinary differential equations. The AM2 model is able to adequately capture the main dynamical behavior of the full anaerobic digestion model ADM1 and has the advantage that a complete analysis for the existence and local stability of its steady states is available. We describe its operating diagram, which is the bifurcation diagram which gives the behavior of the system with respect to the operating parameters represented by the dilution rate and the input concentrations of the substrates. This diagram, is very useful to understand the model from both the mathematical and biological points of view.


Introduction
The anaerobic digestion is a complex process in which organic material is converted into biogas (methane) in an environment without oxygen. Anaerobic digestion enables the water industry to treat waste water as a resource for generating energy and recovering valuable by-products. The complexity of the anaerobic digestion process has motivated the development of complex models, such as the widely used Anaerobic Digestion Model No. 1 (ADM1) [6]. This model has a large number of state variables and parameters. It is impossible to obtain an analytical characterization of the steady states and to describe the operating diagram, that is to say, to identify the asymptotic behaviour of existing steady-states as a function of chemostat operating parameters (substrates inflow concentrations and dilution rate). To the knowledge author, only numerical investigations are available [9].
Due to the analytical intractability of the full ADM1, work has been made towards the construction of simpler models that preserve biological meaning whilst reducing the computational effort required to find mathematical solutions of the model equations, to obtain a better understanding of the anaerobic digestion process. The most simple model of the chemostat with only one biological reaction, where one substrate is consumed by one microorganism is well understood [19,22,29]. However such models are too simple to encapsulate the essence of the anaerobic digestion process.
More realistic models of anaerobic digestion are two-step models, with a cascade of two biological reactions, where one substrate S 1 is consumed by one microorganism X 1 , to produce a product S 2 that serves as the main limiting substrate for a second microorganism X 2 as schematically represented by the following reaction scheme: where µ 1 and µ 2 are the kinetics of the reactions and k i are pseudo-stoichiometric coefficients associated to the bioreactions. An important contribution on the modelling of anaerobic digestion as a two-step is presented by Bernard The two-step models studied in [7,8,28] present a commensalistic relationship between the microorganisms. According to Stephanopolous [30], the commensalism is characterized by the fact that the second population (the commensal population) benefits for its growth from the first population (the host population) while the host population is not affected by the growth of the commensal population and hence, the first population can grow without the second one. On the contrary, the two-step models studied in [11,26,35] present a syntrophic relationship between the microorganisms: the first population is affected by the growth of the second population, and hence no population can grow without the other. For more details and information on commensalism and syntrophy the reader is referred to [10,13,24,25,26,30,31] and the references therein.
Another important and interesting extensions of the two-step anaerobic digestion models are the mathematical models, which include syntrophy and substrate inhibition, considered by Weedermann et al. [33,34] and the three-step models, which consist in introducing an additional microorganism and substrate in a two-step syntrophic model, considered by Wade et al. [32] and Sari and Wade [27].
In this paper we will consider the two-step model of Bernard et al. [8], with general growth functions as in [7,28], denoted here after AM2, and we describe its operating diagram. The operating diagram has the operating parameters as its coordinates and the various regions defined in it correspond to qualitatively different asymptotic behavior. A two-step model has three operating parameters that are the input concentration of substrate for each reaction and the dilution rate. These parameters are control parameters since they are under the control of the experimenter. Apart from these three parameters, that can vary, all other parameters have biological meaning and are fixed depending on the organisms and substrate considered.
Therefore the operating diagram is the bifurcation diagram that shows how the system behaves when we vary the control parameters. As it was claimed by Smith and Waltman in their monography on the mathematical theory of the chemostat (see [29], p. 252), the operating diagram is probably the most useful answer for the discussion of the behavior of the model with respect of the parameters. This diagram shows how robust or how extensive is the parameter region where some asymptotic behavior occur.
This bifurcation diagram is very useful to understand the model from both the mathematical and biological points of view. Its importance for ecological modeling was emphasized by De Freitas et al. [15] and for bioreactors by Pavlou [23]. These authors attributed its introduction to Jost et al. [20], who studied the dynamics of predator and prey interactions in a chemostat. This diagram is often constructed both in the biological literature [15,20,23,28,32,35] and the mathematical literature [1,4,9,11,12,14,19,21,26,27,33,34].
The operating diagram of the AM2 model was only partially described by Sbarciog et al. [28]. In this paper we give a complete description of the diagram. AM2 model can have up to six steady state. Its operating diagram presents nine regions according to the steady state and their stability, that can exist in each region. The operating diagram summarizes the effect of the operating conditions on the long-term dynamics of the AM2 model and shows six type of behavior visualized in the figures by six different colors. Since AM2 model has three operating parameters, and it is not easy to visualize regions in the three-dimensional operating parameter space, two of the operating parameters are used as coordinates of the operating diagram and the effect of the third parameter are shown in a series of operating diagrams.
This paper is organized as follows: in section 2, we present the mathematical model and recall the necessary and sufficient condition of existence and local, and global stability of its steady states. Next, in section 3, we present the operating diagram in the three-dimensional operating parameters space, in sections 4 and 5 we present the operating diagrams in two-dimensional operating parameters space when one of the parameters is kept fixed. In section 6, we present some bifurcations diagram, with the dilution rate as the bifurcation parameter. Then, we conclude by discussing our results in section 7. Proofs and Tables are given in the appendix.

Mathematical model
We consider the AM2 model of anaerobic digestion given in [8], which takes the form of a two-step reactions (1) where, in the first step, the organic substrate S 1 is consumed by the acidogenic bacteria X 1 and produces a substrate S 2 (Volatile Fatty Acids), while, in the second step, the methanogenic population X 2 consumes S 2 and produces methane. Let D be the dilution rate, S 1in and S 2in the concentrations of influent substrate S 1 and S 2 , respectively. The dynamical equations of the model take the form: where k i are pseudo-stoichiometric coefficients associated to the bioreactions and α ∈ [0, 1] is a parameter allowing us to decouple the HRT (Hydraulic Retention Time) and the SRT (Solid Retention Time). In [8], the kinetics µ 1 and µ 2 are of Monod and Haldane type, respectively: The mass flow of the methane production, denoted by Q CH4 , is proportional to the microbial activity, see [5]: where k 4 is the coefficient in (1). In this model the biogas is simply a product of the biological reactions and it has no feedback on the dynamical equations (2).
Following [7,28], we will consider (2) with general C 1 kinetics functions µ 1 and µ 2 satisfying the following qualitative properties: Hypothesis 2. µ 2 (0) = 0, µ 2 (+∞) = 0 and there exists S M 2 > 0 such that µ 2 (S 2 ) > 0 for 0 < S 2 < S M 2 , and µ 2 (S 2 ) < 0 for S 2 > S M 2 . As it is usual in the mathematical theory of the chemostat, see for instance [25], we can use a change of variables that reduces the pseudo-stochiometric coefficients k i to 1. Indeed, the linear change of variables where However, since the stoichiometric coefficients have their own importance for the biologist and since we are interested in giving the biologist a useful tool for the understanding of the role of the operating parameters, we do not make this reduction and we present the results in the original model (2). This model can have at most six steady states, labeled below as in [7]: • E 0 1 , where X 1 = 0 and X 2 = 0: the washout steady state where acidogenic and methanogenic bacteria are extinct.

Operating diagram
Let us consider the surfaces Γ i , i = 1 · · · 6, defined by Table 4. Notice that S 1 * 2 (D) < S 2 * 2 (D) for 0 < D < D 2 and equality holds for D = D 2 . Similarly H 1 (D) < H 2 (D) for 0 < D < min(D 1 , D 2 ), and equality holds for D = min(D 1 , D 2 ). Therefore, the Γ i surfaces separate the operating space (D, S 1in , S 2in ) into nine regions, denoted I k , k = 0 · · · 8, and defined in Table 5. These regions of the operating parameters space (D, S 1in , S 2in ) are corresponding to different system behaviors, as stated in the following result.  Table 5: Definitions of the nine regions corresponding to the nine cases in [7].

Case of [7] Region Definition
Proposition 2. Assume that Hypotheses 1 and 2 hold. The existence and stability properties of the steady states of (2) are given in Table 6, where the regions I k , k = 0 · · · 8 are defined in Table 5.
Proof. The proof is given in Appendix A.2.
Remark 2. In in Figs. 3, 4, 5, 7 and 8 presenting operating diagrams, a region is colored according to the color in Table 6. Each color corresponds to different asymptotic behavior: • Red for the washout of both species, that is, the steady state E 0 1 is Globally asymptotically stable (GAS), which occurs in region I 0 .
• Blue for the washout of acidogenic bacteria while methanogenic bacteria are maintained, that is, the steady state E 1 1 is GAS, which occurs in region I 1 .
• Cyan for the bistability of E 0 1 and E 1 1 which are both (locally) stable. This behavior occurs in region I 1 . Depending on the initial condition the system can go to the washout of both species or the washout of only the acidogenic bacteria.
• Yellow for the washout of methanogenic bacteria while acidogenic bacteria are maintained, that is the steady state E 0 2 is GAS, which occurs in region I 3 .
• Green for the global asymptotic stability of the positive steady state E 1 2 , which occur in I 4 and I 6 . These regions differ only by the existence, in the second region, of the unstable boundary steady state E 1 1 . • Pink for the bistability of E 0 2 and E 1 2 which are both locally asymptotically stable. This behavior occurs in regions I 5 , I 7 and I 8 . These regions differ only by the possible existence of the unstable boundary steady states E 1 1 or E 2 1 . Depending on the initial condition the system can go to the washout of methanogenic bacteria or the coexistence of both species.  Table 11. Compare with Fig. 4 of [28] The operating diagram highly depends on the shape of Γ 4 and Γ 5 surfaces, that is to say, on the behaviors of functions H i , i = 1, 2, defined in Table 1. Notice that these functions are defined on (0, min(D 1 , D 2 )) and H 1 is increasing, since it is the sum of two increasing functions. We have For the limits at right of the domain of definition of these functions, we must distinguish two cases: • When D 1 < D 2 , the functions H i , i = 1, 2 are defined on (0, D 1 ) and Two qualitatively different sub-cases can be distinguished: either H 2 is decreasing on (0, D 2 ) or it is not monotonous. Since H 2 is decreasing near the extremities of its definition interval, a typical example is where it is decreasing, then increasing and then decreasing.
Therefore there are three cases summarized in Table 7 and illustrated in Fig. 1. The role of H i -functions, in the description of the operating diagram, has already been highlighted, see Fig. 4 in [28], where cases D 2 < D 1 and D 1 < D 2 are distinguished.
Since the surfaces Γ i , i = 1 · · · 6, which are the boundaries of the various regions have been derived analytically, the operating diagrams can be drawn qualitatively in each of these cases. Instead of giving a general qualitative description of the operating diagram, and without loss of generality, we present the specific examples shown in Fig. 1. These examples are obtained with the Monod and Haldane functions 3. Notice that these functions satisfy Hypotheses 1 and 2. Therefore, the results of Propositons 1 and 2 apply. The analytical expressions of the auxiliary functions defined in Table 1 and needed in the defintions of the regions I k of the operating diagrams are given in Table 14, in the particular case of functions 3. The biological parameter values used in the figures are given in Table 11. For the sake of practical applicability, these parameter values were chosen in a range that can be found in the literature [7,8].
For the biological parameter values corresponding Fig. 1(a), the surfaces Γ i , i = 1 · · · 6 are shown in Fig. 2. It is difficult to visualize the regions I k , k = 0 · · · 8 of the three-dimensional operating diagram. We can have a better understanding of these regions by showing cuts along 2 dimensional planes where one of the operating parameters is kept constant. For instance, if D is kept constant, we obtain then the operating diagram in the 2-dimensional plane (S 1in , S 2in ). These operating diagrams are described in section 4. If S 2in is kept constant, we obtain then the operating diagram in the 2-dimensional plane (D, S 1in ). These operating diagrams are described in section 5. The intersections of the surfaces Γ i , i = 1 · · · 5 with a plane where D is kept constant are straight lines: vertical line for Γ 1 , horizontal lines for Γ 2 and Γ 3 and oblique lines for Γ 4 and Γ 5 , see Table 12. These straight lines separate the operating parameter plane (S 1in , S 2in ) in up to nine regions I k , k = 0 · · · 8. Since the curves are straight lines, the regions of the operating diagram are very easy to picture. We begin by considering the case where D 2 < D 1 corresponding to Figs. 1(a) and 1(b).

Operating diagram when
The cuts at D constant of the 3-dimensional operating diagram shown in Fig. 2 and corresponding to Fig. 1(a), are shown in Fig. 3. The regions are colored according to the colors in Table 6. For the clarity of the picture all straight lines Γ i are plotted in black. Fig. 3 shows the following features.
For 0 < D < D 2 all regions exist, see Fig. 3(a). For increasing D, the vertical line Γ 1 defined by S 1in = S * 1 (D) moves to the right and tends towards the vertical line defined by S 1in = S * 1 (D 2 ). At the same time, the horizontal lines Γ 2 and Γ 3 , defined by S 2in = S 1 * 2 (D) and S 2in = S 2 * 2 (D), respectively, move towards each other and tend toward the horizontal line defined by S 2in = S M 2 , so that the regions I 1 , I 4 and I 5 shrink an disappear, see Fig. 3 For D 2 < D < D 1 only regions I 0 and I 3 appear, see Figs. 3(c) and 3(d). For increasing D, the vertical line Γ 1 defined by S 1in = S * 1 (D) moves to the right and tends towards infinity. For D ≥ D 1 only region I 0 appears. The cuts D constant of the 3-dimensional operating diagram corresponding to Fig. 1(b), are shown in Fig. 4. This figure has the same qualitative characteristics as Fig. 3: presence of all regions when 0 < D < D 2 as shown in Fig. 4(a); disappearance of all regions except regions I 0 and I 3 , when D = D 2 , as shown in the transition from Fig. 4(b) to Fig. 4(c); disappearance of region I 3 , when D ≥ D 1 , as shown in 4(d). The cuts D constant of the 3-dimensional operating corresponding to Fig. 1(c), are shown in Fig. 5. The regions are colored according to the colors in Table 6. Fig. 5 shows the following features.

Operating diagram when
For 0 < D < D 1 all regions appear, see Fig. 5(a). For increasing D, the vertical line Γ 1 defined by S 1in = S * 1 (D) moves to the right and tends towards infinity. At the same time, the horizontal lines Γ 2 and Γ 3 , defined by S 2in = S 1 * 2 (D) and S 2in = S 2 * 2 (D), respectively, move towards each other, as depicted in Fig. 5(b), and tend towards the horizontal lines defined by S 2in = S 1 * 2 (D 1 ) and S 2in = S 2 * 2 (D 1 ), respectively, as depicted in Fig. 5(c).
For D = D 1 , the operating diagram changes dramatically: all regions I 3 , I 4 and I 5 , I 6 , I 7 and I 8 have disappeared since they are located to the right of the vertical Γ 1 which tends toward infinity, when D tends to D 1 , as depicted in Fig. 5(c).
For D 1 ≤ D < D 2 only regions I 0 , I 1 , and I 2 appear. For increasing D, the horizontal lines Γ 2 and Γ 3 , defined by S 2in = S 1 * 2 (D) and S 2in = S 2 * 2 (D), respectively, move towards each other and tend toward the horizontal line defined by S 2in = S M 2 , so that the regions I 1 shrinks an disappear, see Fig. 5(d). For D = D 2 the operating diagram changes dramatically, since regions I 1 and I 2 disappear and region I 0 invades the whole operating plan. For D ≥ D 2 only region I 0 appears.    Γ k , k = 1 · · · 6, separate the operating parameter plane (D, S 1in ) in up to nine regions I k , k = 0 · · · 8. We begin by considering the case where D 2 < D 1 corresponding to Figs. 1(a) and 1(b).

Operating diagram when
The cuts at S 2in constant of the 3-dimensional operating diagram shown in Fig. 2 and corresponding to Fig. 1(a), are shown in Fig. 6. The regions are colored according to the colors in Table 6. Fig. 6 shows the following features.
For S 2in = 0, only the regions I 0 , I 3 , I 4 and I 5 exist, see Fig. 6(a). For 0 < S 2in < S M 2 , Γ 2 curve appears, giving birth to I 1 , I 6 and I 7 regions, see Fig. 6(b). For increasing S 2in , Γ 4 and Γ 5 curves are translated downwards, while the vertical line Γ 2 moves to the right and tends towards the vertical line Γ 6 , as S 2in tends to S M 2 . For S 2in = S M 2 , Γ 4 curve disappears, while Γ 2 becomes equal to Γ 6 , so that I 4 and I 5 regions have disappeared, see Fig. 6(c). For S 2in > S M 2 , Γ 3 curve appears, giving birth to I 2 and I 8 regions, see Fig. 6(d). For increasing S 2in , the vertical line Γ 3 moves to the left, while Γ 5 curve is translated downwards.
The cuts S 2in constant of the 3-dimensional operating diagram corresponding to Fig. 1(b), are shown in Fig. 7. This figure has the same qualitative characteristics as Fig. 6: presence of only I 0 , I 3 , I 4 and I 5 regions when S 2in = 0, see Fig. 7(a); appearance of I 1 , I 6 and I 7 regions when 0 < S 2in < S M 2 , see Fig. 7(b); disappearance of I 4 and I 5 regions when S 2in = S M 2 , see Fig. 7(c); appearance of I 2 and I 8 regions when S 2in > S M 2 , see Fig. 7(d). It should be noticed that in Fig. 7, the region of global asymptotic stability of the positive steady state E 1 2 (the Green region I 4 ∪ I 6 ) presents the very surprising property that there exists a range of values for the operating parameters S 1in and S 2in such that the system can go from the bistability region (the Pink region I 5 ∪ I 7 ), to the global asymptotic stability region, when the dilution rate D increases. Indeed, the boundary Γ 5 of Green and Pink regions has an increasing part, with respect to parameter D. Therefore, near this part of Γ 5 , as S 1in is kept constant and D increases the system goes from I 5 to I 4 , see Fig. 7(a) and 7(b), or goes from I 7 to I 6 , see Fig. 7(c).
This possibility of globally stabilizing the system, which presents bistability, is surprising since the global stability of the positive steady state is more likely obtained by decreasing D rather than increasing it. This unespecated behavior was first observed in a slightly different two-step model, where the first kinetics is of Contois type [17]. This behavior is investigated in [18].
It is worth-noting that this unexpected behavior can occur only for suitable values of the biological parameters. For instance, in Fig. 6, where all biological parameters are the same as in Fig. 7, excepted that m 1 is changed from m 1 = 0.5 to m 1 = 0.6, the behavior does not occur and a transition from Pink region to Green region is possible only by decreasing D.

Operating diagram when D 1 < D 2
The cuts at S 2in constant of the 3-dimensional operating diagram corresponding to Fig. 1(c), are shown in Fig. 8. The regions are colored according to the colors in Table 6. Since D 1 < D 2 there exists a value S 0 2 < S M 2 such that µ 2 (S 0 2 ) = αD 1 . Fig. 8 shows the following features. For S 2in = 0, only regions I 0 , I 3 , I 4 and I 5 appear, see Fig. 8(a). For 0 < S 2in < S 0 2 , Γ 2 curve appears, giving birth to I 1 , I 6 , I 7 regions, see Fig. 8(b). For increasing S 2in , Γ 4 and Γ 5 curves are translated downwards, while the vertical line Γ 2 moves to the right and tends towards the common vertical asymptote D = D 1 for curves Γ 1 , Γ 4 and Γ 5 , as S 2in tends to S 0 2 . In the limit S 2in = S 0 2 , the very tiny region I 3 (in Yellow on the figure) located between curves Γ 1 and Γ 4 , together with I 4 and I 5 regions have disappeared.
For S 0 2 < S 2in ≤ S M 2 , only regions I 0 , I 1 , I 6 and I 7 exist. For increasing S 2in , the vertical line Γ 2 moves to the right and tends towards Γ 6 as S 2in tends to S M 2 , see Fig. 8(c). For S 2in > S M 2 , Γ 3 curve appears, giving birth to I 2 and I 8 regions, see Fig. 8(d). For increasing S 2in , the vertical line Γ 3 moves to the left while Γ 5 curve is translated downwards.
It should be noticed that as the case (B), it is seen in Fig. 8 that the region of global asymptotic stability of the positive steady state E 1 2 (the Green region I 4 ∪ I 6 ) presents the property that there exists a range of values for the operating parameters S 1in and S 2in such that the system can go from the bistability region (the Pink region I 5 ∪ I 7 ), to the global asymptotic stability region, when the dilution rate D increases.

Bifurcations
The surfaces Γ k , k = 1 · · · 6, are the borders of the regions in the operating parameters space (D, S 1in , S 2in ) on which bifurcations occur, while the steady states change their stability. In codimension-one bifurcations, only transcritical and saddle node bifurcations can be encountered, as stated in the following result. Proposition 3. The bifurcations of the steady states of (2) arising on the boundaries of regions I k , k = 0 · · · 8, are listed in Table 8.
Proof. The proof is given in Appendix A.3. Table 8: Codimension-one bifurcations along subsets of surfaces Γ k and the corresponding cases in [7]: Transcritical bifurcations (TB) and Saddle Node bifurcations (SNB) occur.
Remark 3. The last column of Table 8 shows the corresponding cases with non hyperbolic steady states given in Theorem 1 of [7]. The case labeled 2.10 in this theorem, where E 0 1 = E 1 1 and E 0 2 = E 2 2 , does not appear in Table  8, since it is a codimension-two bifurcation arising along Γ 2 ∩ Γ 5 . The bifurcations along Γ 1 , corresponding to the condition S 1in = S * 1 (D) were not analyzed in [7]. In Theorem 1 of [7] only the cases S 1in < S * 1 (D) and S 1in > S * 1 (D) were considered.
To have a better understanding of the nature of the bifurcations of steady states, let us consider the dilution rate D as the bifurcation parameter. Throughout this section, we assume that biological parameters are fixed as in Fig. 7(a), corresponding to case (b) of Fig. 1 and S 2in = 0. We now fix the operating parameter S 1in at various typical values, as depicted in the horizontal lines shown in Fig. 9, and plot one-parameter bifurcation diagrams in D, with X i , i = 1, 2, on the y-axis, see Fig. 10  Recall that the curve Γ 5 separating the Pink and Green regions is the curve of the function S 1in = k1 k2 H 2 (D). Case (B) corresponds to a function H 2 which is decreasing, then increasing, then decreasing. For the considered biological parameters values, the function H 2 (D) attains its minimum for D min 0.72 and its maximum for D max = 0.81 and satisfies We fix three typical values S 1in = 13, S 1in = 13.3 and S 1in = 14, corresponding to the three horizontal lines shown in Fig. 9. We begin with the case where S 1in = 14. Since S 1in > 13.37, as it is seen in Fig. 9, with increasing D, there is a transition from I 4 to I 5 for D = D 3 ≈ 0.5917, then from I 5 to I 3 for D = D 2 0.8186, then from I 3 to I 0 for D = D 4 0.8696. The bifurcation values D 2 , D 3 and D 4 are defined by The bifurcation value D 4 corresponds to a transcritical bifurcation of E 0 2 and E 0 1 ; D 2 corresponds to a saddle node bifurcation of E 1 2 and E 2 2 and D 3 corresponds to a transcritical bifurcation of E 0 2 and E 2 2 . The plot of X 1 and X 2 components of all existing steady states with respect of D is shown in Fig. 10. Solid lines and dotted lines correspond to stable and unstable steady states respectively. Since S 2in the steady states E 1 1 and E 2 1 cannot exist. On Fig. 10(a), for Consider now the case where S 1in = 13. This case corresponds to the surprising situation where we can go from the bistability region (colored in Pink) to the global asymptotic stability region (colored in Green), when the dilution rate D increases. Since 12.57 < S 1in < 13.11, as it is seen in Fig. 9, with increasing D, there is a transition from I 4 to I 5 for D = D 7 ≈ 0.6526, then from I 5 to I 4 for D = D 6 0.7844, then from I 4 to I 3 for D = D 5 0.8184, then from The bifurcation value D 4 corresponds to a transcritical bifurcation of E 0 2 and E 0 1 ; D 5 corresponds to a transcritical bifurcation of E 1 2 and E 0 2 and D 6 and D 7 correspond to transcritical bifurcations of E 0 2 and E 2 2 . The plot of X 1 and X 2 components of all existing steady states with respect of D is shown in Fig. 11. Solid lines and dotted lines correspond to stable and unstable steady states respectively. On Fig. 11(a), for 0 < D < D 7 and D 6 < D < D 5 the X 1 -component of E i 2 , i = 0, 2, is colored in Red, with Green dots, showing the stability of E 1 2 and the instability of E 0 2 . For D 7 < D < D 6 , the X 1 -component of E i 2 , i = 0, 1, 2, is colored in Red and Green, with Blue dots, showing the bistability of E 0 2 and E 1 2 and the instability of E 1 2 . For D 5 < D < D 4 , the X 1 -component of E 0 2 is colored in Green, showing the stability of E 0 2 . On Fig. 11(b) and 11(c), for 0 < D < D 7 and D 6 < D < D 5 , the X 2 = 0-component of E 0 j , j = 1, 2, is colored with Green and Black dots, showing the instability of E 0 2 and E 0 1 ; For D 7 < D < D 6 and D 5 < D < D 4 it is colored in Green, with Black dots, showing the stability of E 0 2 and the instability of E 0 1 . For D > D 4 it is colored in Black showing the stability of E 0 1 . Consider now the case where S 1in = 13.3. This case corresponds also to the situation where we can go from the bistability region (colored in Pink) to the global asymptotic stability region (colored in Green), when the dilution rate D increases. Since 13.11 < S 1in < 13.37, as it is seen in Fig. 9, with increasing D, there is a transition from I 4 to I 5 for D = D 10 ≈ 0.6304, then from I 5 to I 4 for D = D 9 0.8050, then from I 4 to I 5 for D = D 8 0.8173, then from The bifurcation value D 4 corresponds to a transcritical bifurcation of E 0 2 and E 0 1 ; D 2 corresponds to a saddle node bifurcation of E 1 2 and E 2 2 and D 8 , D 9 and D 10 correspond to transcritical bifurcations of E 0 2 and E 2 2 . The plot of the X 2 component of all existing steady states with respect of D is shown in Fig. 12. Solid lines and dotted lines correspond to stable and unstable steady states respectively. Since two magnifications are necessary to represent all bifurcations, the plot of the X 1 component is omitted in Fig. 12

Discussion
Our main contribution is to present the operating diagram and to show how it depends on the biological parameters.
The parameter space of model (2), where µ 1 and µ 2 are given by (3) is twelve dimensional: nine biological or physical parameters (m 1 , m 2 , K 1 , K 2 ,K I , k 1 , k 2 , k 2 and α) and three operating parameters (D, S 1in and S 1in ). The former parameters are called biological parameters since they depend on the organisms, and substrate considered. These parameters are measurable in the laboratory. In contrast, the later parameters are called operating parameters since they are under the control of the experimenter.
Exploring all of the twelve dimensional parameter space is almost possible. Fixing the biological parameters and constructing the operating diagram is a powerful answer for the discussion of the behavior of the model with respect of the parameters. Therefore our approach to handle the question of the dependence with respect of the parameters of the model is to split the question in two intermediary questions. First we fix the biological parameters and present the operating diagram. Second we explore how the operating diagram varies when the biological parameters are changed. For instance, Figs. 6, Fig. 7 and Fig. 8 show how the operating diagram changes when the biological parameter m 1 is changed.
The operating diagrams shown in the figures summarize the effect of the operating conditions on the long-term dynamics of the AM2 model and shows six type of behavior: 1) the washout of the two populations (regions colored in Red), 2) the washout of the first population while the second population is maintained (regions colored in Blue), 2) the occurrence of these two behaviors, according to initial conditions (regions colored in Cyan), 4) the washout of the second population while the first is maintained (regions colored in Yellow), 5) the persistence of both populations (regions colored in Green) and finally 6) the occurrence of these two behaviors according to initial conditions (regions colored in Pink).
In the operating diagrams shown in Figs. 6(a), 7(a) and 8(a), obtained for S 2in = 0, only regions I 0 , I 3 , I 4 and I 5 exist, that is to say, the steady states E i 1 , i = 1, 2 without acidogenic bacteria, cannot exist. This property is in accordance with the fact that the system being commensalistic, and without input concentration S 2in , it is impossible for the commensal population (the methanogenic bacteria) to survive if the host population (the acidogenic bacteria) is washed out.
The operating diagram shows how robust or how extensive is the parameter region where coexistence occurs, where the corresponding steady state is GAS, where the steady states, with extinction both or one of the population, is stable and where it is unstable.

A Proofs
A.1 Proof of Proposition 1 Table 9: Necessary and sufficient conditions of existence and local stability of the steady states of (2) obtained in [7]. S * 1 (D), S i * 2 (D) and S * 2in (D, S 1in , S 2in ) are defined in Table 1.

Existence conditions Stability conditions
The proof follows from [7]. It is seen from Proposition 1 of [7] that the steady states are given by Table 2, where S * 1 , S i * 2 , S * 2in , X * 1 , X i 2 and X i * 2 are defined in Table 1. The necessary and sufficient conditions of existence of the steady state given in Proposition 1 of [7] are summarized in the second column of Table 9. The necessary and sufficient conditions of local stability of these steady states, obtained in Table A.1 of [7], are summarized in the third column of Table 9.
Let us prove the following result which shows that the existence conditions of X i * 2 , i = 1, 2 steady states given in Table  9, can be stated using the functions H i (D), i = 1, 2, defined in Table 1. These functions were considered also by [28]. Lemma 1. The conditions S * 2in (D, S 1in , S 2in ) = S i * 2 (D) and S * 2in (D, S 1in , S 2in ) < S i * 2 (D), for i = 1, 2, are equivalent to the conditions S 2in + k2 k1 S 1in = H i (D) and S 2in + k2 k1 S 1in < H i (D), for i = 1, 2, respectively.
Proof. The result follows from the definitions of S * 2in (D, S 1in , S 2in ) and H i (D), given in Table 1.
That is to say S 2in + k2 k1 S 1in = H i (D). The proof for the inequality is the same.
Therefore, the results in Table 9 are equivalent to those in Table 3 which completes the proof of Proposition 1.

A.2 Proof of Proposition 2
The proof follows from [7]. The existence and stability conditions of the steady states of (2) given in Table 9 depend only on the relative positions of the values of S 1in and S * 1 (D) and of the values of S 1 * 2 (D), S 2 * 2 (D), S 2in , and S * 2in (D, S 1in , S 2in ). Actually, as stated in Theorem 1 of [7], we can distinguish nine cases, according to the relative positions of these numbers. These cases are summarized in Table 10.
The cases 1.1, 1.2 and 1.3 correspond to the regions I 0 , I 1 and I 2 respectively, defined in Table 5. Now we use Lemma 1 to show that the remaining six cases 2.1 to 2.6 correspond to the six regions I 3 to I 8 defined in Table 5. Table 10: The 9 cases of existence and stability of steady states of (2) obtianed in [7], where S and U stand for stable and unstable respectively.
Since S 2in < S * 2in the case 2.1 corresponds to the condition S * 2in < S 1 * 2 which is equivalent, using Lemma 1, to S 2in + k2 k1 S 1in < H 1 (D). Therefore the case 2.1 corresponds to the region I 3 defined in Table 5. Using again Lemma 1, the condition S 1 * 2 < S * 2in < S 2 * 2 in the case 2.2 is equivalent to H 1 (D) < S 2in + k2 k1 S 1in < H 2 (D) and the condition S * 2in > S 2 * 2 in the case 2.3 is equivalent to S 2in + k2 k1 S 1in > H 2 (D). Therefore the cases 2.2 and 2.3 correspond to the regions I 4 and I 5 respectively, defined in Table 5. Using similar arguments we show that the cases 2.4, 2.5 and 2.6 correspond to the regions I 6 , I 7 and I 8 respectively, defined in Table 5.
Excepted for cases 1.3, 2.3, 2.5 and 2.6 of bistability, the system (2) has a unique globally asymptotically stable (GAS) steady state. Therefore, in the case 1.1, E 0 1 is GAS; in the case 1.2, E 1 1 is GAS, in the case 2.1, E 0 2 is GAS, and in the cases 2.2 and 2.4, E 1 2 is GAS. In the case 1.3, E 2 1 is a saddle point whose attractive manifold is a 3-dimensional hyper-surface surface which separates the phase space of (2) into the basins of attractions of the stable steady states E 0 1 and E 1 1 . In the cases 2.3, 2.5 and 2.6, E 2 2 is a saddle point whose stable manifold is a 3-dimensional hyper-surface which separates the phase space of (2) into the basins of attractions of the stable steady states E 0 2 and E 1 2 . For details and complements on the global behaviour, see section 2.4 of [7]. This completes the proof of Proposition 2.

B Tables
In this section, we give several tables that are used in the paper. In Table 11, we provide the biological parameter values used in the figures. Tables 12 and 13, we give the description of the intersection of the Γ k surfaces with a two dimesnional operating plane where D or S 2in is kept constant respectively. In Table 14, we present the functions defined in Table 1 in the particular case of the Monod and Haldane growth function 3.  [7] and corresponding to the figures.
The whole plane if D = D 2 Empty if D = D 2 Table 13: The intersections of the Γ k surfaces, k = 0 · · · 8 with a (D, S 1in ) plane, where S 2in is kept constant.