Tewfik Sari
email: tewfik.sari@irstea.fr

Jérôme Harmand
email: harmand@inra.supagro.fr

A model of a syntrophic relationship between two microbial species in a chemostat including maintenance

Keywords: Microbial ecosystems, Syntrophic relationship, Maintenance, Stability, Food chains 2010 MSC: 34D20, 92D25, 92C45

Many microbial ecosystems can be seen as microbial 'food chains' where the different reaction steps can be seen as such: the waste products of the organisms at a given reaction step are consumed by organisms at the next reaction step.

In the present paper we study a model of a two-step biological reaction with feedback inhibition, which was recently presented as a reduced and simplified version of the anaerobic digestion model ADM1 of the International Water Association (IWA). It is known that in the absence of maintenance (or decay) the microbial 'food chain' is stable. In a previous study, using a purely numerical approach and ADM1 consensus parameter values, it was shown that the model remains stable when decay terms are added. However, the authors could not prove in full generality that it remains true for other parameter values. In this paper we prove that introducing decay in the model preserves stability whatever its parameters values are and for a wide range of kinetics.

Introduction

Two-step models are commonly used to describe microbial systems, which take the form of a cascade of two biological reactions where one substrate S 0 is consumed by one microorganism X 0 to produce a product S 1 that serves as the main limiting substrate for a second microorganism X 1 as schematically represented by the following reaction scheme:

S 0 f0 -→ X 0 + S 1 , S 1 f1 -→ X 1
where f 0 and f 1 are the growth functions that may depend on several substrates.

The substrate and biomass concentrations in this two-step model evolve according to the four-dimensional dynamical system of ODEs

                     dS 0 dt = D S in 0 -S 0 -f 0 (•) X 0 dX 0 dt = -αDX 0 + Y 0 f 0 (•) X 0 -a 0 X 0 dS 1 dt = D S in 1 -S 1 + Y 2 f 0 (•) X 0 -f 1 (•) X 1 dX 1 dt = -αDX 1 + Y 1 f 1 (•) X 1 -a 1 X 1 (1)
where Y 0 are the Yield coefficients. Substrate S 0 and S 1 are introduced with an input concentration S in 0 and S in 1 respectively, and at dilution rate D. Depending on the technology used to confine the reactions, the coefficient α ≤ 1 is not necessarily equal to 1 and 1 -α represents the proportion of biomass which is retained in the reactor. This model includes the maintenance (or decay) terms a 0 and a 1 . Maintenance, in its most general assertion, is the consumption of energy for all processes other than growth: it is modelled either by adding a negative term on the substrate dynamic without associating it to growth or by considering a decay term on the biomass dynamics, as in [START_REF] Alcaraz-Gonzalez | Software sensors for highly uncertain WWTPs : a new approach based on interval observers[END_REF]. For more information about the modelling of maintenance, the reader is referred to [START_REF] Ni | Model-based characterization of endogenous maintenance, cell death and predation processes of activated sludge in sequencing batch reactors[END_REF]. These models present the advantage of being complex enough to capture important process properties while being simple enough to be mathematically studied.

When the growth function f 0 depends only on the substrate S 0 and the growth function f 1 depends only on the substrate S 1 , that is

f 0 (•) = f 0 (S 0 ), f 1 (•) = f 1 (S 1 ), ( 2
)
the system is known as commensalistic: one species grows on the product of another one [START_REF] Reilly | Stability of commensalistic systems[END_REF][START_REF] Stephanopoulos | The dynamics of commensalism[END_REF]. The system has a cascade structure: solve the first and second equations for S 0 , X 0 , and then use this result is the remaining equations to find S 1 , X 1 . Consequently S 0 and X 0 are the same in pure and mixed culture experiments. The number of steady-states and their stability as a function of model inputs and parameters may be investigated [START_REF] Benyahia | Bifurcation and stability analysis of a two step model for monitoring anaerobic digestion processes[END_REF][START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF][START_REF] Sbarciog | Determination of appropriate operating strategies for anaerobic digestion systems[END_REF].

When f 0 depends on both substrates S 0 and S 1 and f 1 depends only on S 1 , that is

f 0 (•) = f 0 (S 0 , S 1 ), f 1 (•) = f 1 (S 1 ) (3)
the system is known as syntrophic. For instance if the first organism is inhibited by high concentrations of the product S 1 , the extent to which the substrate S 0 is degraded by the organism X 0 depends on the efficiency of the removal of the product S 1 by the bacteria X 1 . The mathematical analysis of such model is more delicate than commensalistic models, see for instance [START_REF] Burchard | Substrate degradation by a mutualistic association of two species in the chemostat[END_REF][START_REF] Kreikenbohm | A mathematical model of syntrophic cocultures in the chemostat[END_REF][START_REF] Kreikenbohm | Bistability in the Chemostat[END_REF][START_REF] Wilkinson | Interactions in a Mixed Bacterial Population Growing on Methane in Continuous Culture[END_REF] and the more recent papers [START_REF] Hajji | A mathematical study of a syntrophic relationship of a model of anaerobic digestion process[END_REF][START_REF] Harvey | Quantifying the effects of the division of labor in metabolic pathways[END_REF][START_REF] Sari | The mathematical analysis of a syntrophic relationship between two microbial species in a chemostat[END_REF][START_REF] Volcke | Steady-state multiplicity of two-step biological conversion systems with general kinetics[END_REF][START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF].

A model of a two-tiered microbial 'food chain' with feedback inhibition, which encapsulates the essence of the anaerobic digestion process was recently proposed [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF]. Anaerobic digestion is a biological process that converts organic matter into a gaseous mixture composed mainly of methane and carbon dioxide through the action of a complex bacterial and archaeal ecosystem. It is often used for the treatment of concentrated wastewaters or to convert the excess sludge produced in wastewater treatment plants into more stable products [START_REF] Ramirez | Modelling microbial diversity in anaerobic digestion through an extended adm1 model[END_REF][START_REF] Van Lier | Anaerobic wastewater treatment[END_REF]. One of its advantages is that the methane produced can be used profitably as a source of energy. It is usually considered that a number of metabolic groups of microorganisms are involved sequentially in several serial and parallel conversion steps to finally produce methane and carbon dioxide. The Anaerobic Digestion Model No. 1 (ADM1) of the IWA Task Group for Mathematical Modelling of Anaerobic Digestion Processes [START_REF] Batstone | The Iwa Anaerobic Digestion Model No 1 (ADM1)[END_REF][START_REF]IWA Task Group for Mathematical Modelling of Anaerobic Digestion Processes[END_REF] is too complex to permit mathematical analysis of its nonlinear dynamics and only numerical investigations are available [START_REF] Bornhöft | Steadystate analysis of the Anaerobic Digestion Model No. 1 (ADM1)[END_REF].

The model of Xu et al. [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF] includes maintenance terms and considers the syntrophic associations between propionate degraders and methanogens. It was shown that the non-trivial steady-state is not necessarily stable. In addition simulation results with the ADM1 consensus values indicate that the positive steady-state is always stable whenever it exists. For the operators of anaerobic wastewater treatment systems the results of Xu et al. [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF] show that the syntrophic associations between propionate degraders and methanogens are inherently stable under realistic environmental conditions. However, the possibility of an unstable positive steady-state was not excluded for other parameter values and the title of [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF], Maintenance affects the stability of a two-tiered microbial 'food chain' ? left unanswered the question of the effects of maintenance from a more general viewpoint. In the present paper, we show that for any values of the parameters the positive steady-state is stable as long as it exists, that is to say, maintenance does not affect the stability of the considered two-tiered microbial 'food chain', see [START_REF] Sari | Maintenance does not affect the stability of a two-tiered microbial 'food chain[END_REF].

The paper is organized as follows. In Section 2, we review the different twostep models that have been proposed in the literature and we recall what useful informations were obtained for applications from their qualitative mathematical analysis. In Section 3 we present the hypothesis on [START_REF] Alcaraz-Gonzalez | Software sensors for highly uncertain WWTPs : a new approach based on interval observers[END_REF][START_REF] Batstone | The Iwa Anaerobic Digestion Model No 1 (ADM1)[END_REF]. In Section 4 we give the description of the steady-state and their stability. In Section 5 we describe the operating diagram. In Section 6 we apply our results to the model of Xu et al. [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF] and we give an answer to open questions on the stability of the positive steady-state of their model. Concluding remarks are given in Section 7. The technical proofs of the results are given in the Appendix.

Commensalism, mutualism and syntrophy

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 1 or a coefficient α in the dilution rate of the biomass is considered or not. In most cases, the models used are not generic in the sense either model parameters are fixed or the growth functions are predefined (Monod, Haldane, etc).

Following Stephanopoulos [START_REF] Stephanopoulos | The dynamics of commensalism[END_REF] we say that 'Two populations of microorganisms which grow in a mixed culture and interact in such a way that one population (the commensal population) depends for its growth on the other population and thus benefits from the interaction while the other population (the host) is not affected by the growth of the commensal population constitutes an example of commensalism'. Reilly [START_REF] Reilly | Stability of commensalistic systems[END_REF] was the first to propose a mathematical study of a pure commensalistic model (1,2) holds and a 0 = a 1 = 0, α = 1. He was interested in explaining surprising oscillations observed within the course of an experiment realized in making Saccharomyces carlsbergensis growing on fructose produced by Acetobacter suboxyduns from mannitol. In particular, he established theoretical conditions involving a feedback from the yeast to the bacteria. In this study, explicit growth functions modelling the proposed feedback were used.

An important contribution on the modelling of anaerobic digestion as a commensalistic system is the model by Bernard et al. [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF]. The authors considered a Monod function for f 0 and a Haldane function for f 1 . Sbarciog et al. [START_REF] Sbarciog | Determination of appropriate operating strategies for anaerobic digestion systems[END_REF] studied this model for α = 1 while the interesting case where 0 < α < 1 and where growth functions were characterized by qualitative properties was studied by Benyahia et al. [START_REF] Benyahia | Bifurcation and stability analysis of a two step model for monitoring anaerobic digestion processes[END_REF]. Prior to these investigations, and regarding the potential of anaerobic systems to produce renewable energy, the study of these models were particularly important for optimizing anaerobic digestion, notably through the synthesis of state observers and control feedback laws (cf. for instance [START_REF] Alcaraz-Gonzalez | Software sensors for highly uncertain WWTPs : a new approach based on interval observers[END_REF][START_REF] Alcaraz-Gonzalez | Application of a Robust Interval Observer to an Anaerobic Digestion Process[END_REF]).

Another fundamental ecological interactions which can be modelled by two-step reaction models with two microorganisms are mutualism and syntrophy.

Mutualism is defined as a situation where two organisms cooperate typically in producing mutually the substrate necessary to the growth of the other [START_REF] Hajji | Association between competition and obligate mutualism in a chemostat[END_REF].

A syntrophic relationship between two organisms refers to growth functions of the form (3) where the species exhibit mutualism but where, in contrast to what happens in a purely symbiotic relationship, one of the species can grow without the other. Important results of these studies were conditions under which a stable coexistence may occur. Wilkinson et al. [START_REF] Wilkinson | Interactions in a Mixed Bacterial Population Growing on Methane in Continuous Culture[END_REF] considered the case of growth functions of the form

f 0 (S 0 , S 1 ) = m 0 S 0 K 0 + S 0 1 1 + S 1 /L 1 , f 1 (S 1 ) = m 1 S 1 K 1 + S 1 (4)
Kreikenbohm et al. [START_REF] Kreikenbohm | A mathematical model of syntrophic cocultures in the chemostat[END_REF] considered the case where f 1 is a Monod function and the growth function f 0 takes the form

f 0 (S 0 , S 1 ) =    m0(S0-S1/L) K0+S0+K1S1 if S 0 -S1/L > 0 0 otherwise (5)
In this case the first organism is unable to grow unless the quotient S 1 /S 0 is small enough, say, S 1 < LS 0 . Burchard [START_REF] Burchard | Substrate degradation by a mutualistic association of two species in the chemostat[END_REF] extended the results of [START_REF] Kreikenbohm | A mathematical model of syntrophic cocultures in the chemostat[END_REF][START_REF] Wilkinson | Interactions in a Mixed Bacterial Population Growing on Methane in Continuous Culture[END_REF] to a large class of more generic growth functions, including the special cases (4) and

(5). He highlighted conditions under which there is persistence or extinction.

El Hajji et al. [START_REF] Hajji | A mathematical study of a syntrophic relationship of a model of anaerobic digestion process[END_REF], motivated by the analysis of the main studied steps of the anaerobic digestion where H 2 -producing acetogens are associated to H 2 -utilizing bacteria, considered the general case where the growth functions (3), satisfy the following properties:

∂f 0 ∂S 0 > 0, ∂f 0 ∂S 1 < 0, df 1 dS 1 > 0 (6)
Another extension was considered by Kreikenbohm et al. [START_REF] Kreikenbohm | Bistability in the Chemostat[END_REF], which considered the case where S 0 appears also in f 1 (•):

f 1 (S 0 , S 1 ) = m 1 S 1 K 1 + S 1 1 1 + S 0 /L 0
The mathematical analysis of this model showed the occurrence of bistability that cannot be observed when f 1 (•) depends only on S 0 . Sari et al. [START_REF] Sari | The mathematical analysis of a syntrophic relationship between two microbial species in a chemostat[END_REF] considered the general situation of a growth function f 1 (•) = f 1 (S 0 , S 1 ), which is increasing in S 1 and decreasing in S 0 and showed, in contrast with the case where f 1 (•) = f 1 (S 1 ) depends only on S 1 , that a multiplicity of positive equilibria can occur. This work was motivated by the study of the influence of the presence of an input term into the dynamics of S 1 again and by the consideration of more general forms for growth rate functions to investigate the association of H 2 -producing acetogens and H 2 -utilizing bacteria. Other models for which

f 0 (•) = f 0 (S 0 , S 1 ) and f 1 (•) = f 1 (S 0 , S 1 )
, exhibiting the multiplicity of positive equilibria can be found in [START_REF] Volcke | Steady-state multiplicity of two-step biological conversion systems with general kinetics[END_REF].

All these studies do not include maintenance terms. This short review of the existing literature shows that under conditions like (6) and without maintenance terms (a 0 = a 1 = 0), the positive steady-state is unique and stable, if its exists [START_REF] Burchard | Substrate degradation by a mutualistic association of two species in the chemostat[END_REF][START_REF] Hajji | A mathematical study of a syntrophic relationship of a model of anaerobic digestion process[END_REF][START_REF] Kreikenbohm | A mathematical model of syntrophic cocultures in the chemostat[END_REF][START_REF] Wilkinson | Interactions in a Mixed Bacterial Population Growing on Methane in Continuous Culture[END_REF]. On the other hand as soon as f 1 (•) = f 1 (S 0 , S 1 ) may depend on S 0 then instability of the positive steady-state can occur [START_REF] Kreikenbohm | Bistability in the Chemostat[END_REF][START_REF] Sari | The mathematical analysis of a syntrophic relationship between two microbial species in a chemostat[END_REF][START_REF] Volcke | Steady-state multiplicity of two-step biological conversion systems with general kinetics[END_REF].

To the best of our knowledge, Xu et al. [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF] were the first to consider the effects of maintenance terms in [START_REF] Alcaraz-Gonzalez | Software sensors for highly uncertain WWTPs : a new approach based on interval observers[END_REF][START_REF] Batstone | The Iwa Anaerobic Digestion Model No 1 (ADM1)[END_REF], in the particular case of the growth functions (4), and S in 1 = 0, α = 1. As mentioned in the introduction these authors were not able to show that the positive steady-state is stable if it exists. In the present paper we will consider the general case [START_REF] Alcaraz-Gonzalez | Software sensors for highly uncertain WWTPs : a new approach based on interval observers[END_REF][START_REF] Batstone | The Iwa Anaerobic Digestion Model No 1 (ADM1)[END_REF] where growth functions satisfy [START_REF] Boer | Food chain dynamics in the chemostat[END_REF] and with maintenance terms (a 0 > 0, a 1 > 0) and S in 1 = 0, α = 1. We will prove that the positive steady-state is stable whenever it exists. Therefore, in this paper we generalize [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF] by allowing a larger class of growth functions, we generalize [START_REF] Wilkinson | Interactions in a Mixed Bacterial Population Growing on Methane in Continuous Culture[END_REF] by allowing a larger class of growth functions and maintenance terms, and we generalize [START_REF] Hajji | A mathematical study of a syntrophic relationship of a model of anaerobic digestion process[END_REF] by allowing maintenance terms. For the applications our results show that that the syntrophic associations between propionate degraders and methanogens are inherently stable for a wide range of kinetics and whatever the parameters values are, not only for the kinetics [START_REF] Benyahia | Bifurcation and stability analysis of a two step model for monitoring anaerobic digestion processes[END_REF] and with the ADM1 consensus values of parameters as shown in [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF].

An important and interesting extension should be mentioned here: Weedermann et al. [START_REF] Weedermann | Mathematical Model of Anaerobic Digestion in a Chemostat: Effects of Syntrophy and Inhibition[END_REF] proposed an 8-dimensional mathematical model, which includes syntrophy and inhibition, both mechanisms considered by Bernard et al. [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF] and by El Hajji et al. [START_REF] Hajji | A mathematical study of a syntrophic relationship of a model of anaerobic digestion process[END_REF]. The effects of maintenance terms are considered by Weedermann et al. [START_REF] Weedermann | Optimal biogas production in a model for anaerobic digestion[END_REF].

The model

In this paper, we study the model (1,3) with α = 1 and S in 1 = 0. We obtain the following system

                     dS 0 dt = D S in 0 -S 0 -f 0 (S 0 , S 1 ) X 0 dX 0 dt = -DX 0 + Y 0 f 0 (S 0 , S 1 ) X 0 -a 0 X 0 dS 1 dt = -DS 1 + Y 2 f 0 (S 0 , S 1 ) X 0 -f 1 (S 1 ) X 1 dX 1 dt = -DX 1 + Y 1 f 1 (S 1 ) X 1 -a 1 X 1 (7)
Notice that we do not assume any specific analytical expression for the growth and inhibition functions. Our analysis will use only the following general assumptions for the growth functions f 0 (S 0 , S 1 ) and f 1 (S 1 ) :

A1 For all S 0 > 0 and S 1 ≥ 0, f 0 (S 0 , S 1 ) > 0 and f 0 (0, S 1 ) = 0.

A2 For all S 1 > 0, f 1 (S 1 ) > 0 and f 1 (0) = 0.

A3 For all S 0 > 0 and S 1 > 0, ∂f 0 ∂S 0 (S 0 , S 1 ) > 0 and ∂f 0 ∂S 1 (S 0 , S 1 ) < 0.

A4 For all S 1 > 0, df 1 dS 1 (S 1 ) > 0.

Hypothesis A1 signifies that no growth can take place for species X 0 without the substrate S 0 . Hypothesis A1 means that the intermediate product S 1 is necessary for the growth of species X 1 . Hypothesis A3 means that the growth rate of species X 0 increases with the substrate S 0 but it is self-inhibited by the intermediate product S 1 . Hypothesis A4 means that the growth of species X 1 increases with intermediate product S 1 produced by species X 0 . Note that this defines a syntrophic relationship between the two species.

To ease the mathematical analysis of the system, we can rescale system [START_REF] Bornhöft | Steadystate analysis of the Anaerobic Digestion Model No. 1 (ADM1)[END_REF] using the following change of variables adapted from [START_REF] Sari | The mathematical analysis of a syntrophic relationship between two microbial species in a chemostat[END_REF]:

s 0 = Y 2 S 0 , x 0 = Y 2 Y 0 X 0 , s 1 = S 1 , x 1 = 1 Y 1 X 1 ,
We obtain the following system

                     ds 0 dt = D(s in 0 -s 0 ) -µ 0 (s 0 , s 1 )x 0 dx 0 dt = -Dx 0 + µ 0 (s 0 , s 1 )x 0 -a 0 x 0 ds 1 dt = -Ds 1 + µ 0 (s 0 , s 1 )x 0 -µ 1 (s 1 )x 1 dx 1 dt = -Dx 1 + µ 1 (s 1 )x 1 -a 1 x 1 (8)
where s in 0 = Y 2 S in 0 and µ 0 and µ 1 are defined by

µ 0 (s 0 , s 1 ) = Y 0 f 0 1 Y 2 s 0 , s 1 and µ 1 (s 2 ) = Y 1 f 1 (s 1 ) (9)
The functions µ 0 and µ 1 are general functions with their own properties. Since the functions f 0 and f 1 satisfy hypotheses A1-A4, it follows from (9) that functions µ 0 and µ 1 satisfy: H1 For all s 0 > 0 and s 1 ≥ 0, µ 0 (s 0 , s 1 ) > 0 and µ 0 (0, s 1 ) = 0.

H2 For all s 1 > 0, µ 1 (S 1 ) > 0 and µ 1 (0) = 0.

H3 For all s 0 > 0 and s 1 > 0, ∂µ 0 ∂s 0 (s 0 , s 1 ) > 0 and ∂µ 0 ∂s 1 (s 0 , s 1 ) < 0.

H4 For all s 1 > 0, dµ 1 ds 1 (s 1 ) > 0.

It should be noticed that (8) was studied in [START_REF] Hajji | A mathematical study of a syntrophic relationship of a model of anaerobic digestion process[END_REF][START_REF] Sari | The mathematical analysis of a syntrophic relationship between two microbial species in a chemostat[END_REF] in the case where maintenance effects are not taken into account, i.e. a 0 = a 1 = 0. We can easily prove that that for every non-negative initial condition, the solution of (8) has non-negative components and is positively bounded and thus is defined for every positive t.

Steady-state and stability analysis

A steady-state of ( 8) is a solution of the following nonlinear algebraic system obtained from (8) by setting the right-hand sides equal to zero:

D(s in 0 -s 0 ) -µ 0 (s 0 , s 1 )x 0 = 0 ( 10
)
-Dx 0 + µ 0 (s 0 , s 1 )x 0 -a 0 x 0 = 0 (11) -Ds 1 + µ 0 (s 0 , s 1 )x 0 -µ 1 (s 1 )x 1 = 0 (12) -Dx 1 + µ 1 (s 1 )x 1 -a 1 x 1 = 0 (13)
A steady-state exists (or is said to be 'meaningful' [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF]) if and only if all its components are non-negative. From equation [START_REF] Harvey | Quantifying the effects of the division of labor in metabolic pathways[END_REF] we deduce that:

x 0 = 0 or µ 0 (s 0 , s 1 ) = D + a 0 ( 14
)
and from equation ( 13) we deduce that:

x 1 = 0 or µ 1 (s 1 ) = D + a 1 (15)
The case x 0 = 0 and x 1 > 0 is excluded. Indeed, as a consequence of (15), we have µ 1 (s 1 ) = D +a 1 and, as a consequence of ( 12), we have Ds 1 +(D +a 1 )x 1 = 0, which is impossible since s 1 ≥ 0 and x 1 > 0. Therefore, three cases must be distinguished: SS0: x 0 = 0, x 1 = 0 where both species are washed out.

SS1: x 0 > 0, x 1 = 0, where species x 1 is washed out while x 0 survives. SS2: x 0 > 0, x 1 > 0, where both species survive.

For the description of the steady-states and their stability, we need the following notations. Since the function s 1 → µ 1 (s 1 ) is increasing, it has an inverse function y → M 1 (y), so that, for all s 1 ≥ 0 and y ∈ [0, sup µ 1 (•))

s 1 = M 1 (y) ⇐⇒ y = µ 1 (s 1 ) (16)
Let s 1 be fixed. Since the function s 0 → µ 0 (s 0 , s 1 ) is increasing, it has an inverse function y → M 0 (y, s 1 ), so that, for all s 0 , s 1 ≥ 0, and y ∈ [0, sup µ 0 (•, s 1 ))

s 0 = M 0 (y, s 1 ) ⇐⇒ y = µ 0 (s 0 , s 1 ) (17)
The inverse functions s 1 = M 1 (y) and s 0 = M 0 (y, s 1 ) can be calculated explicitly in the case of the Monod growth functions [START_REF] Stephanopoulos | The dynamics of commensalism[END_REF] considered in Section 6, see formulas (24,(25)). We define the functions:

F 0 (D) = M 0 (D + a 0 , 0) F 1 (D) = M 1 (D + a 1 ) + M 0 (D + a 0 , M 1 (D + a 1 )) (18)
Steady-state Existence condition Stability condition SS0 Always exists Notice that F 1 (D) > F 0 (D) for all D ≥ 0, as long as they are both defined with the exception F 1 (0) = F 0 (0), which holds if and only if a 0 = a 1 = 0. Now, we can describe the steady-states of (8).

s in 0 < F 0 (D) SS1 s in 0 > F 0 (D) s in 0 < F 1 (D) SS2 s in 0 > F 1 (D) Always Stable
Proposition 1. Assume that assumptions H1-H4 hold. Then (8) has at most three steady-states:

• SS0= s 0 = s in 0 , x 0 = 0, s 1 = 0, x 1 = 0 It always exists. It is stable if and only if s in 0 < F 0 (D).
• SS1= s 0 , x 0 = D D+a0 s in 0 -s 0 , s 1 = s in 0 -s 0 , x 1 = 0 where s 0 is the solution of equation µ 0 (s 0 , s in 0 -s 0 ) = D + a. It exists if and only if s in 0 > F 0 (D). It is stable if and only if s in 0 < F 1 (D).

• SS2= s 0 , x 0 = D D+a0 s in 0 -s 0 , s 1 , x 1 = D D+a1 s in 0 -s 0 -s 1 where s 1 = M 1 (D + a 1 ) and s 0 = M 0 (D + a 0 , M 1 (D + a 1 )). It exists if and only if s in 0 > F 1 (D). It is stable if it exists.
The proof is given in the Appendix.

Notice that SS1 exists as soon as SS0 becomes unstable and SS2 exists as soon as SS1 becomes unstable. One concludes that for any value of the operating parameters, there is always one, and only one, steady-state which is stable. The results are summarized in Table 1. When decay effects are not taken into account, i.e. a 0 = a 1 = 0, the system can be reduced to a planar system and global stability results can be obtained [START_REF] Hajji | A mathematical study of a syntrophic relationship of a model of anaerobic digestion process[END_REF][START_REF] Sari | The mathematical analysis of a syntrophic relationship between two microbial species in a chemostat[END_REF]: for any pair of operating parameters, there is always one, and only one, steady-state which is globally asymptotically stable.

(a

) (b) s in 0 s in 0 D D J 1 J 0 J 2 Γ 0 Γ 1 u0 u1 J 2 J 1 J 0 Γ 0 Γ 1

Operating diagram

The operating diagram shows how the system behaves when we vary the two control parameters S in 0 and D. Let F 0 (D) and F 1 (D) be the functions defined by [START_REF] Reilly | Stability of commensalistic systems[END_REF]. The curve Γ 0 of equation s in 0 = F 0 (D) is the border which makes SS0 unstable and at the same time SS1 exists (the red curve on Fig. 1). The curve Γ 1 of equation s in 0 = F 1 (D) is the border which makes SS1 unstable and at the same time SS2 exists (the blue curve on Fig. 1).

The curves Γ 0 and Γ 1 separate the operating plane (s in 0 , D) in three regions, as shown in Fig. 1, labelled J 0 , J 1 and J 2 . The results of Prop. 1 are summarized in Table 2 which shows the existence and stability of the steady-states SS0, SS1 and SS2 in the regions J 0 , J 1 and J 2 of the operating diagram.

The values u 0 and u 1 plotted on the figure are obtained as follows:

u 0 = F 0 (0) = M 0 (a 0 , 0), u 1 = F 1 (0) = β + M 0 (a 0 , β) , β = M 1 (a 1 ) (19)
If a 0 ≥ sup s0>0 µ 0 (s 0 , 0), F 0 (0) is not defined and we let u 0 = +∞. In this case the regions J 1 and J 2 are empty. If a 1 < sup s1>0 µ 1 (s 1 ) or a 0 ≥ sup s0>0 µ 0 (s 0β, β), F 1 (0) is not defined and we let u 1 = +∞. In this case the region J 2 is empty. When maintenance effects are not taken into consideration, then u 0 = u 1 = 0 and we have

F 0 (D) = M 0 (D, 0), F 1 (D) = M 1 (D) + M 0 (D, M 1 (D)) Condition Region SS0 SS1 SS2 s in 0 < F 0 (D) (s in 0 , D) ∈ J 0 S F 0 (D) < s in 0 < F 1 (D) (s in 0 , D) ∈ J 1 U S F 1 (D) < s in 0 (s in 0 , D) ∈ J 2 U U S

A two-tiered microbial 'food chain'

The model considered in [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF] involves a two-tiered microbial 'food chain' with feedback inhibition, consisting of a propionate degrader and a hydrogenotrophic methanogen. The propionate degrader produces hydrogen which inhibits its own growth. Using the notations of ADM1 the model can be written as

                     dS pro dt = D (S pro,in -S pro ) -f 0 (S pro , S H2 ) X pro dX pro dt = -DX pro + Y pro f 0 (S pro , S H2 ) X pro -k dec,pro X pro dS H2 dt = -DS H2 + 0.43 (1 -Y pro ) f 0 (S pro , S H2 ) X pro -f 1 (S H2 ) X H2 dX H2 dt = -DX H2 + Y H2 f 1 (S H2 ) X H2 -k dec,H2 X H2 (20)
where S pro and X pro are propionate substrate and biomass concentrations; S H2 and X H2 are those for hydrogen; Y pro and Y H2 are the Yield coefficients and 0.43 (1 -Y pro ) represents the part which goes to hydrogen substrate. Both growth functions take Monod form with an hydrogen inhibition for the first one

f 0 (S pro , S H2 ) = k m,pro S pro K s,pro + S pro 1 1 + S H 2 K I,H 2 , f 1 (S H2 ) = k m,H2 S H2 K s,H2 + S H2 (21)
Here, apart from the two operating (or control) parameters, which are the inflowing propionate concentration S pro,in and the dilution rate D, that can vary, all others have biological meaning and are fixed depending on the organisms and substrate considered [see 30, Table 1]. The aim of Xu et al. The basic results of the analysis of [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF] are: for any pair of values of operating parameters, at most one steady-state is stable. When one of the decay terms is not taken into account, i.e. k dec,pro = 0 or k dec,H2 = 0 in [START_REF] Sari | Maintenance does not affect the stability of a two-tiered microbial 'food chain[END_REF], there is always one and only one steady-state which is stable and SS2 is stable as long as it exists. When both decay effects are present, i.e. k dec,pro > 0 and k dec,H2 > 0 in [START_REF] Sari | Maintenance does not affect the stability of a two-tiered microbial 'food chain[END_REF], the authors were not able to check all the Routh-Hurwitz criteria for SS2. They claimed that SS2 is not necessarily stable in theory when it exists and they established numerically that with the ADM1 parameters values, SS2 is stable as long as it exists. However they did not give any values for the biological parameters for which, under some operating parameters, SS2 becomes unstable. As a consequence of Proposition 1, we can say that, for all values of the parameters, SS2 is stable whenever it exits, which actually gives an answer to the questions asked by [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF] in their paper.

More precisely, using the following simplified notations in ( 20)

S 0 = S pro , S in 0 = S pro,in , S 1 = S H2 , X 0 = X pro , X 1 = X H2 Y 0 = Y pro , Y 1 = Y H2 , Y 2 = 0.43 (1 -Y pro ) , a 0 = k dec,pro , a 1 = k dec,H2
and using the rescaling ( 9) and the biological parameters in ( 21) we obtain the model ( 8) with the following growth function:

µ 0 (s 0 , s 1 ) = m 0 s 0 K 0 + s 0 1 1 + s 1 /K i , µ 1 (s 1 ) = m 1 s 1 K 1 + s 1 ( 22
)
where

m 0 = Y 0 k m,pro , K 0 = Y 2 K s,pro , K i = K I,H2 m 1 = Y 1 k m,H2 , K 1 = K s,H2 (23)
Let us describe our results in the particular case [START_REF] Burchard | Substrate degradation by a mutualistic association of two species in the chemostat[END_REF][START_REF] Stephanopoulos | The dynamics of commensalism[END_REF]. Notice that the growth functions [START_REF] Stephanopoulos | The dynamics of commensalism[END_REF] satisfy Assumptions H1-H4, so that Proposition 1 holds. In this case the inverse functions M 1 (y) and y → M 0 (y, s 1 ) of the functions µ 1 (s 1 ) and s 0 → µ 0 (s 0 , s 1 ) can be calculated explicitly: we have

y ∈ [0, m 1 ) → M 1 (y) = K 1 y m 1 -y , ( 24
)
y ∈ 0, m 1 1 + s 1 /K i → M 0 (y, s 1 ) = K 0 y m0 1+s1/Ki -y (25)
Therefore, the functions F 1 (D) and F 2 (D) defined by ( 18) are given explicitly by

F 0 (D) = K 0 (D + a 0 ) m 1 -D -a 0 F 1 (D) = K 1 (D + a 1 ) m 1 -D -a 1 + K 0 (D + a 0 ) m 0 1 + K1(D+a1) (m1-D-a1)Ki -D -a 0 (26)
Notice that F 0 is defined on [0, m 1 -a 0 ) and F 1 is defined on [0, D + ) with D + < m 1 -a 0 . On the other hand, the solution s 0 of equation µ 0 (s 0 , s in 0 -s 0 ) = D+a 0 , which is used in SS1, is simply the positive solution of the quadratic equation:

m 0 s 0 = (D + a 0 )(K 0 + s 0 ) 1 + s in 0 -s 0 K i ( 27
)
As a corollary of Proposition 1 we have the following result.

Proposition 2. Assume that µ 0 and µ 1 are given by [START_REF] Stephanopoulos | The dynamics of commensalism[END_REF]. Let F 0 (D) and F 1 (D) be defined by [START_REF] Volcke | Steady-state multiplicity of two-step biological conversion systems with general kinetics[END_REF]. Then (8) has at most three steady-states

• SS0= s 0 = s in 0 , x 0 = 0, s 1 = 0, x 1 = 0 It always exists. It is stable if and only if s in 0 < F 0 (D). • SS1= s 0 , x 0 = D D+a0 s in 0 -s 0 , s 1 = s in 0 -s 0 , x 1 =
0 where s 0 is the positive solution of the quadratic equation [START_REF] Weedermann | Mathematical Model of Anaerobic Digestion in a Chemostat: Effects of Syntrophy and Inhibition[END_REF]. It exists if and only if s in 0 > F 0 (D). If it exists then it is stable if and only if

s in 0 < F 1 (D). • SS2= s 0 , x 0 = D D+a0 s in 0 -s 0 , s 1 , x 1 = D D+a1 s in 0 -s 0 -s 1 where s 1 = K 1 (D + a 1 ) m 1 -D -a 1 , s 0 = K 0 (D + a 0 ) m 0 1 + s1 Ki -D -a 0 It exists if and only if s in 0 > F 1 (D). It is stable if it exists.
As a consequence of this result we obtain the results of [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF]. To make the comparison possible the reader is advised on the main difference between our approach and [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF]: we use the rescaling ( 9) and hence work with the growth functions [START_REF] Stephanopoulos | The dynamics of commensalism[END_REF], while Xu et al. For the numerical simulations we will use the nominal values of Table 3 obtained from Table 1 of [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF] by using the formulas [START_REF] Van Lier | Anaerobic wastewater treatment[END_REF] and a 0 = k dec,pro ,

m 0 d -1 0.52 K 0 kg COD/m 3 0.124 m 1 d -1 2.10 K 1 kg COD/m 3 2.5 10 -5 K i kg COD/m 3 3.5 10 -6 a 0 d -1 0.02 a 1 d -1 0.02
a 1 = k dec,H2
. For these values of the parameters, the values u 0 and u 1 are very small, see Fig. 2. Notice that the scaling on the two coordinates in Fig. 2 are different from those of Fig. 2 of [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF], since these authors used another rescaling.

Discussion

Following [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF], we considered a two-tiered 'food chain' with feedback inhibition, which is a generalized model describing the syntrophic interaction of a propionate degrader and a hydrogenotropic methanogen. In the absence of (a)

(b) s in 0 s in 0 D D J 1 J 0 J 2 Γ 0 Γ 1 J 2 J 1 J 0 Γ 0 Γ 1 u0 u1
Figure 2: Operating diagram of the model ( 8)-( 22). (a) The model was parametrised with the ADM1 consensus values listed in Table 3. (b) A magnification showing the values u 0 = 49.6 10 -4 , u 1 = 53.1 10 -4 defined by [START_REF] Sari | The mathematical analysis of a syntrophic relationship between two microbial species in a chemostat[END_REF].

maintenance these authors proved that this two-tiered 'food chain' is always stable. When maintenance is included in the model they were not able to check the Routh-Hurwitz criteria, and since the possibility of having at least one pair of complex eigenvalues with positive real parts is not theoretically excluded, they concluded that Hopf bifurcation can originate from SS2 [see 30, Appendix B].

However, using the consensus parameters of ADM1 and numerical simulations, they have shown that the model of the methanogenic two-tiered propionatehydrogen food chain is always stable [see 30, Section 6.2]. In this work we have generalized the model of the two-tiered 'food chain' of [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF] by considering generic growth functions and we established the stability of the generalized model with maintenance terms.

In [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF], the authors point out that introducing decay or maintenance in the classical predator-prey models results in instability and chaos [START_REF] Kooi | Chaotic behaviour of a predator-prey system in the chemostat[END_REF]. For more details on food-chains in the chemostat the reader may consult [START_REF] Boer | Food chain dynamics in the chemostat[END_REF][START_REF] Van Voorn | Ecological consequences of global bifurcations in some food chain models[END_REF][START_REF] Vayenas | Chaotic dynamics of a food web in a chemostat[END_REF].

Therefore, they observed that, in spite of the fundamental differences between their 'food chain' and the classical predator-prey models, the same intrinsic effect of maintenance on the stability of the food chain is observed [see 30,

Section 7]: When maintenance is included in its description, the two-tiered generalized 'food chain' is not necessarily stable in theory. The results obtained in the present paper indicate that the two-tiered generalized 'food chain' is always stable, so that it is fundamentally different from the classical predatorprey model.

It should be noticed that the rescaling used by [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF] gives a dimensionless model. However, our present rescaling [START_REF] Hajji | Association between competition and obligate mutualism in a chemostat[END_REF] does not yield a dimensionless model.

The new variables s 0 , x 0 , s 1 and x 1 have the same dimensions as the original variables S 0 , X 0 , S 1 and X 1 . The original growth functions [START_REF] Sbarciog | Determination of appropriate operating strategies for anaerobic digestion systems[END_REF] are transformed by our rescaling [START_REF] Hajji | Association between competition and obligate mutualism in a chemostat[END_REF] in the growth functions [START_REF] Stephanopoulos | The dynamics of commensalism[END_REF] with the same scale imbalance in the half-saturation rates, see Table 3. We cannot benefit from the dimensionless rescaling used by Xu et al. [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF], because this rescaling uses some kinetics parameters of the specific growth functions (21) while we work with general unspecified growth functions. The benefit of our rescaling ( 9) is that it permits to fix the constant yields parameters Y 0 , Y 1 and Y 2 in (7) to 1, as shown by the rescaled model [START_REF] Burchard | Substrate degradation by a mutualistic association of two species in the chemostat[END_REF].

We were successful in checking the Routh-Hurwitz criteria because we work with general growth functions (defined by their qualitative properties given in assumptions A1-A4) and our computations are not encumbered by the specific form of the growth functions considered by Xu et al. [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF]. These authors noticed [see 30, Section 7] that direct application of symbolic analysis programs, such as Maple or Mathematica, did not provide adequate solutions for the stability of the system. Actually we used the symbolic analysis program Maple to verify that the coefficients β i in the expression of the term

f 1 f 2 f 3 -f 2 1 f 4 -f 2 3
given in Appendix D are correct. It should be noticed that [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF] have claimed [see 30,

Remark 1] that their method is still effective for other growth functions. Our main contribution was to believe them and to try to solve the problem with general growth functions.

In the model [START_REF] Bornhöft | Steadystate analysis of the Anaerobic Digestion Model No. 1 (ADM1)[END_REF] considered in this work, the first species X 0 uses the substrate S 0 for its growth and produces a substrate S 1 consumed by the second species X 1 for its growth. The substrate S 1 produced by the first species inhibits its own growth, that is, the growth function f 0 (S 0 , S 1 ) is decreasing with respect to S 1 . In practice, and in many complex models as the ADM1, it hap-pens that the second species is also inhibited by the first substrate. Thus, it is interesting to consider the case where the second species is inhibited by the substrate S 0 , namely that f 1 (S 0 , S 1 ) also depends on S 0 and is decreasing with respect to S 0 . It has been shown by Sari et al. [19] that the introduction of this last inhibiting relationship in the model completely changes the model properties while maintenance was not considered. In particular, the modified model exhibits multiplicity of positive steady-states. However, it should be stressed that these results were very general: whether this instability occurs for realistic environmental conditions or not is under investigation.

Another interesting question, which is the object of a future work, is to consider an input term S in 1 in (7), as well as a coefficient α < 1 in the dilution rate of the biomass, as it was the case in the general setting of (1). For instance if S in 1 > 0 then there exists an additional steady-state where X 0 = 0 is washed out and X 1 > 0 does not go to extinction.

Appendix A. Stability analysis

We give the proof of Prop. 1. A steady-state (s 0 , x 0 , s 1 , x 1 ) of ( 8) is a solution of the set of algebraic equations [START_REF] Hajji | A mathematical study of a syntrophic relationship of a model of anaerobic digestion process[END_REF][START_REF] Harvey | Quantifying the effects of the division of labor in metabolic pathways[END_REF][START_REF]IWA Task Group for Mathematical Modelling of Anaerobic Digestion Processes[END_REF][START_REF] Kooi | Chaotic behaviour of a predator-prey system in the chemostat[END_REF]. The local stability of each steady-state depends on the sign of the real parts of the eigenvalues of the corresponding Jacobian matrix for the system (8). This is the matrix of the partial derivatives of the right hand side with respect to the state variables evaluated at the given steady-state (s 0 , x 0 , s 1 , x 1 ), that is:

J =         -D -Ex 0 -µ 0 F x 0 0 Ex 0 µ 0 -D -a 0 -F x 0 0 Ex 0 µ 0 -D -F x 0 -Gx 1 -µ 1 0 0 Gx 1 µ 1 -D -a 1         (A.1)
where

E = ∂µ 0 ∂s 0 (s 0 , s 1 ) > 0, F = - ∂µ 0 ∂s 1 (s 0 , s 1 ) > 0, G = dµ 1 ds 1 (s 1 ) > 0
The eigenvalues of J are the roots of its characteristic polynomial det(J -λI).

Notice that we have used the opposite sign for the partial derivative F = SS1 exists if and only if this equation has a solution in the interval (0, s in 0 ). The function s 0 → ψ(s 0 ) = µ 0 (s 0 , s in 0 -s 0 ) is increasing since its derivative dψ ds0 = ∂µ0 ∂s0 -∂µ0 ∂s1 > is positive. Using ψ(0) = 0 and ψ(s in 0 ) = µ(s in 0 , 0) we conclude that equation (C.1) has a solution in the interval (0, s in 0 ) if and only if ψ(s in 0 ) = µ(s in 0 , 0) > D + a 0 , that is to say condition (B.1) holds. The condition of existence of SS1 is then equivalent to the condition of instability of SS0. Evaluated at SS1, the Jacobian matrix (A.1) becomes:

J =         -D -Ex 0 -D -a 0 F x 0 0 Ex 0 0 -F x 0 0 Ex 0 D + a 0 -D -F x 0 -µ 1 0 0 0 µ 1 -D -a 1         Its characteristic polynomial is: det(J-λI) = (λ-µ 1 +D+a 1 )(λ+D) λ 2 + [D + (E + F )x 0 ] λ + (D + a 0 )(E + F )x 0
Its eigenvalues are λ 1 = µ 1 -D -a 1 , λ 2 = -D and λ 3 and λ 4 are the roots of the following quadratic equation:

λ 2 + [D + (E + F )x 0 ] λ + (D + a 0 )(E + F )x 0 = 0 Since λ 3 λ 4 = (D + a 0 )(E + F )x 0 > 0 and λ 3 + λ 4 = -[D + (E + F )x 0 ] < 0,
the real parts of λ 3 and λ 4 are negative. So for being stable it must be λ 1 < 0.

Therefore SS1 is stable if and only if

µ 1 (s in 0 -s 0 ) < D + a 1 , where s 0 is the solution of (C.1) (C.2)
Since the function s 1 → µ 1 (s 1 ) is increasing, we have the following equivalence

µ 1 (s in 0 -s 0 ) < D + a 1 ⇐⇒ s 0 < s in 0 -M 1 (D + a 1 )
Since the function s 0 → ψ(s 0 ) = µ 0 s 0 , s in 0 -s 0 is decreasing, we deduce that ψ (s 0 ) > ψ s in 0 -M 1 (D + a 1 ) . Since s 0 be the solution of (C. Using [START_REF] Ramirez | Modelling microbial diversity in anaerobic digestion through an extended adm1 model[END_REF] we obtain s 0 = M 0 (D + a 0 , M 1 (D + a 1 )). As a result of ( 10) and ( 12)

x 0 = D D + a 0 s in 0 -s 0 , x 1 = D D + a 1 s in 0 -s 0 -s 1
SS2 exists if and only if s in 0 > s 0 + s 1 , that is β 4 = 4(Hx 0 + Gx 1 ) 2 + 2a 0 Hx 0 + 2a 1 Gx 1 β 3 = 2(Hx 0 + Gx 1 ) 3 + 4EG(Hx 0 + Gx 1 )x 0 x 1 +5a 0 H 2 x 2 0 + (a 0 + a 1 )(3E + 5F )Gx 0 x 1 + 5a 1 G 2 x 2 1 β 2 = 4EG(Hx 0 + Gx 1 ) 2 x 0 x 1 +3a 0 H 3 x 3 0 + (a 0 E + 2a 1 E + 6a 0 H + 3a 1 F )GHx 2 0 x 1 +(2a 0 E + a 1 E + 3a 0 F + 6a 1 H)G 2 x 0 x 2 1 + 3a 1 G 3 x 3 1 +a 2 0 F (F + 2E)x 2 0 + (a 0 Ex 0 -a 1 Gx 1 ) 2 + 2a 0 a 1 GF x 0 x 1 β 1 = 2E 2 G 2 (Hx 0 + Gx 1 )x 2 0 x 2 1 + (4a 0 + a 1 )EGH 2 x 3 0 x 1 +(a 0 + a 1 )(3E + 5F )EG 2 x 2 0 x 2 1 + (a 0 + 4a 1 )EG 3 x 0 x 3 1 +a 2 0 (3E 2 + 3EF + F 2 )F x 3 0 + a 0 (2a 0 E + a 0 F + 2a 1 F )GF x 2 0 x 1 +(Ex 0 + Gx 1 )(a 0 Ex 0 -a 1 Gx 1 ) 2 + (2a 0 a 1 + a 2 1 )G 2 F x 0 x 2 1 β 0 = (a 0 + a 1 )E 2 G 2 (Hx 0 + Gx 1 )x 2 0 x 2 1 + a 2 0 (2E + F )EF Gx Therefore, SS2 is stable as long as it exists.

Figure 1 :

 1 Figure 1: Operating diagram without (a) and with (b) maintenance effects. The values u 0and u 1 are defined by[START_REF] Sari | The mathematical analysis of a syntrophic relationship between two microbial species in a chemostat[END_REF]

[ 30 ]

 30 was to study the stability of the steady-states of the model[START_REF] Sari | Maintenance does not affect the stability of a two-tiered microbial 'food chain[END_REF][START_REF] Sbarciog | Determination of appropriate operating strategies for anaerobic digestion systems[END_REF] while varying the two operating (or control) parameters D and S pro,in . The system[START_REF] Sari | Maintenance does not affect the stability of a two-tiered microbial 'food chain[END_REF][START_REF] Sbarciog | Determination of appropriate operating strategies for anaerobic digestion systems[END_REF] can have at most three steady-states: a trivial solution where both populations are washed out (SS0), a solution where X H2 is washed out while X pro survives (SS1) and a positive solution where both populations survive (SS2). The local stability of each steady-state was tested by linearisation around the steady-state values of the variables.

  [START_REF] Xu | Maintenance affects the stability of a two-tiered microbial chain' ?[END_REF] use a dimensionless rescaling. Despite this difference, both approaches are equivalent and hence must give the same results. Our quadratic equation (27) used in the description of SS1 is the same as their quadratic equation (A.1), or the quadratic equation without numbering preceding equation (B.1) [see 30, Appendix A and B]. Parameters Units Nominal Value

1 Using ( 16 )

 116 1), ψ (s 0 ) = µ 0 s 0 , s in 0 -s 0 = D + a 0 Therefore, the condition (C.2) of stability of SS1 is equivalent to:D + a 0 < µ 0 s in 0 -M 1 (D + a 1 ), M 1 (D + a 1 ) (C.3) Since the function s 0 → µ 0 (s 0 , M 1 (D + a 1 )) is increasing, and using (17), the condition (C.3) is equivalent tos in 0 -M 1 (D + a 1 ) < M 0 (D + a 0 , M 1 (D + a 1 ))which is, according to (18), equivalent tos in 0 < M 1 (D + a 1 ) + M 0 (D + a 0 , M 1 (D + a 1 )) =: F 1 (D) Appendix D. SS2At SS2, x 0 = 0, x 1 = 0. As a consequence of (14) and (15) s 0 and s 1 are solutions of the set of equations µ 0 (s 0 , s 1 ) = D + a 0 , µ 1 (s 1 ) = D + a we obtain s 1 = M 1 (D + a 1 ) and s 0 is a solution of equation µ 0 (s 0 , M 1 (D + a 1 )) = D + a 0 (D.1)

s in 0 > M 1 ( 4 where f 1 =f 2 = 2 f 3 =f 1 f 2 -f 3 = 2D 3 + α 2 D 2 + α 1 D + α 0 where α 2 = 4 ( 1 Thus f 1 f 2 -f 3 On the other hand we have f 1 f 2 f 3 -f 2 1 f 4 -f 2 3 = β 5 D 5 + β 4 D 4 + β 3 D 3 + β 2 D 2 + β 1 D + β 0 where β 5 = 2 (

 14122323322124123355443322152 D + a 1 ) + M 0 (D + a 0 , M 1 (D + a 1 )) =: F 1 (D) Evaluated at SS2, the Jacobian matrix (A.1) becomes:Ex 0 D + a 0 -D -F x 0 -Gx 1 -D -a 1Its characteristic polynomial is:det(J -λI) = λ 4 + f 1 λ 3 + f 2 λ 2 + f 3 λ + f Gx 1 + (E + F )x 0 + 2D EGx 0 x 1 + (2D + a 0 )(E + F )x 0 + (2D + a 1 )Gx 1 + D (2D + a 0 + a 1 )EGx 0 x 1 + D(D + a 0 )(E + F )x 0 + D(D + a 1 )Gx 1 f 4 = (D + a 0 )(D + a 1 )EGx 0 x 1 Hence f i > 0 for i = 1 • • • 4 (D.2)Since the quantity E+F occurs so often in the computations, we use the notationH = E + F . Straightforward calculations show that: Hx 0 + Gx 1 ) α 1 = 2(Hx 0 + Gx 1 ) 2 + a 0 Hx 0 + a 1 Gx 1 α 0 = EG(Hx 0 + Gx 1 )x 0 x 1 + a 0 H 2 x 2 0 + (a 0 + a 1 )F Gx 0 x 1 + a 1 G 2 x 2 Hx 0 + Gx 1 )

3 0 x 1 +(a 2 0 + a 2 1 )EF G 2 x 2 0 x 2 1 + (a 0 1 Thus f 1 f 2 f 3 -f 2 1 f 4 -f 2 3 > 0 (D. 4 )

 121101404 Ex 0 -a 1 Gx 1 ) 2 EGx 0 xAccording to (D.2), (D.3) and (D.4) the Routh-Hurwitz criteria are satisfied.

Table 1 :

 1 Existence and local stability of steady-states.

Table 2 :

 2 Existence and local stability of steady-states. The letter S (resp. U) means stable (resp. unstable). No letter means that the steady-state does not exist.

Table 3 :

 3 Nominal parameters values.

Acknowledgements

This work was supported by the Treasure network (https://project.inria.fr/treasure/).

The authors warmly thank

Matthew Wade for his comments. This work was finalized during the stimulating workshop, Microbial ecology and mathematical modelling, organized in December 2014 at the Centre Interfacultaire Bernoulli (CIB) of the Ecole Polytechnique Fédérale de Lausanne (EPFL) in Switzerland, with funding of the Swiss National Science Foundation.

Appendix B. SS0

At SS0, x 0 = 0, x 1 = 0. As a result of [START_REF] Hajji | A mathematical study of a syntrophic relationship of a model of anaerobic digestion process[END_REF] and [START_REF]IWA Task Group for Mathematical Modelling of Anaerobic Digestion Processes[END_REF], s 0 = s in 0 and s 1 = 0. SS0 always exists. Evaluated at SS0, the Jacobian matrix (A.1) becomes

For being stable we need λ 1 < 0. Therefore SS0 is unstable if and only if

Since the function s 0 → µ 0 (s 0 , 0) is increasing, and using [START_REF] Ramirez | Modelling microbial diversity in anaerobic digestion through an extended adm1 model[END_REF] we have the following equivalence

Therefore, according to ( 18), (B.1) is equivalent to s in 0 > F 0 (D).

Appendix C. SS1

At SS1, x 0 = 0, x 1 = 0. As a consequence of ( 14) µ 0 (s 0 , s 1 ) = D + a 0 . As a result of ( 10) and ( 12) D(s in 0 -s 0 ) = µ 0 (s 0 , s 1 )x 0 and Ds 1 = µ 0 (s 0 , s 1 )x 0 Hence x 0 = D D+a0 s in 0 -s 0 and D(s in 0 -s 0 ) = Ds 1 , so that s 0 + s 1 = s in 0 . Therefore s 0 is a solution of equation µ 0 (s 0 , s in 0 -s 0 ) = D + a 0 (C.1)


Julian Per

Benoit Becker
email: benoit.celse@ifp.fr

Celse

Denis Guillaume

Victor Costa

Luc Bertier

Emman- Uelle Guillon

Gerhard Pirngruber

A Continuous Lumping Model for Hydrocracking on a Zeolite Catalysts: Model Development and Parameter Identification

Keywords: Hydrocracking, Kinetic Modeling, Continuous Lumping, Sensitivity Study

Hydrocracking of non-conventional and non-crude feedstocks is becoming an increasingly active area of research. This study is focused on the development of a hydrocracking model based on traditional VGO residues.

. There are two main motivations for developing hydrocracking models: 1) obtaining robust predictions of the performance of HCK units, which can be used to guide process design and optimization, and 2) obtaining a more nuanced understanding of the underlying physico-chemical processes to guide basic research and development of new catalysts and/or process designs.

Accurate and robust simulation of residue (i.e. >370°C cut) conversion, yield structure in terms of standard cuts, and

Introduction

The rapid decline in conventional light crude oil resources [START_REF] Ward | Hydrocracking processes and catalysts[END_REF], combined with the increasing demand for middle distillates (i.e. Kerosene and Gasoil cuts) [START_REF] Motaghi | Slurry-phase hydrocracking-possible solution to refining margins: Opportunity crudes require more hydrogen addition to upgrade orphan product streams into higher-value "clean" products[END_REF][START_REF] Butler | Maximize liquid yield from extra heavy oil: Next-generation hydrocracking processes increase conversion of residues[END_REF] have made Hydrocracking an increasingly attractive option for upgrading of Vacuum Gas Oil (VGO) residue. The VGO cut is characterized by a high True Boiling Point (i.e. high percentage of >370°C cut) [START_REF] Riazi | Characterization and Properties of Pretroleum Fractions[END_REF]. It frequently contains relatively high amounts of organic nitrogen and sulfur, which are undesirable in the final products (ie.g. diesel, fuel oil). The purpose of hydrocracking units is the transformation of VGO to valuable middle distillate fractions, conforming to ever more stringent product quality specifications [START_REF] Rana | A review of recent advances on process technologies for upgrading of heavy oils and residua[END_REF].

The hydrocracking (HCK) process involves the breaking up of large hydrocarbons by β-scission, hydrogenation of aromatic rings, as well as PCP-isomerization. The removal of organic nitrogen and sulfur by hydrodenitrogenation (HDN) and hydrodesulfurization (HDS) reactions are another important function of a HCK unit [START_REF] Rana | A review of recent advances on process technologies for upgrading of heavy oils and residua[END_REF]. A hydrocracking unit often uses a zeolite catalyst with an acid and a metallic site [START_REF]Catalyse acido-basique[END_REF]. A good general description of the mechanism is given by Weitkamp [START_REF] Weitkamp | Catalytic Hydrocracking-Mechanisms and Versatility of the Process[END_REF].

The reaction mechanism of such bifunctional catalysts has been the object of a number of past [START_REF] Scherzer | Hydrocracking Science and Technology[END_REF][START_REF] Martens | Attempts to Rationalize the Distribution of Hydrocracking Products. I Qualitative Description of the Primary Hydrocracking Modes of Long Chain Paraffins in open[END_REF] and recent studies [START_REF] Henry | Vacuum gas oil hydrocracking performance of bifunctional Mo/Y zeolite catalysts in a semi-batch reactor[END_REF][START_REF] Burnens | Hydrocracking reaction pathways of 2,6,10,14tetramethylpentadecane model molecule on bifunctional silica-alumina and ultrastable Y zeolite catalysts[END_REF][START_REF] Roussel | Mechanisms of n-decane hydrocracking on a sulfided NiW on silica-alumina catalyst[END_REF][START_REF] Francis | Design of improved hydrocracking catalysts by increasing the proximity between acid and metallic sites[END_REF][START_REF] Pellegrini | Modeling of Fischer-Tropsch products hydrocracking[END_REF][START_REF] Browarzik | Hydrocracking process of n-alkanes by continuous kinetics[END_REF][START_REF] Botchwey | Product Selectivity during Hydrotreating and Mild Hydrocracking of Bitumen-Derived Gas Oil[END_REF]. Many of these studies are focused on model feeds, such as Fisher-Tropsch wax [START_REF] Pellegrini | Modeling of Fischer-Tropsch products hydrocracking[END_REF] or n-alkane mixtures [START_REF] Browarzik | Hydrocracking process of n-alkanes by continuous kinetics[END_REF]. Maya residue [START_REF] Sanchez | Kinetic Model for Moderate Hydrocracking of Heavy Oils[END_REF] and Athabasca bitumen derived VGO [START_REF] Botchwey | Product Selectivity during Hydrotreating and Mild Hydrocracking of Bitumen-Derived Gas Oil[END_REF] are among the industrial feedstocks used for model development.. Hydrocracking is considered to be an extremely versatile process [START_REF] Ward | Hydrocracking processes and catalysts[END_REF][START_REF] Weitkamp | Catalytic Hydrocracking-Mechanisms and Versatility of the Process[END_REF]. A large number of different feedstocks can be converted into a broad range of products by carefully choosing the catalyst and adjusting operating conditions [START_REF] Ward | Hydrocracking processes and catalysts[END_REF][START_REF] Dry | High quality diesel via the Fischer-Tropsch process -a review[END_REF][START_REF] Choudhary | Hydrocracking: A Review[END_REF][START_REF] Calemma | Middle distillates from hydrocracking of FT waxes: Composition, characteristics and emission properties[END_REF]. simulated distillation (SIMDIS) is required by a hydrocracking model. The hydrogen consumption, generation of gases (C3 & C4), as well as a number of product qualities should also be well predicted. An advanced HCK model should be validated for a wide range of operating conditions and for different feeds in order to have the explanatory capabilities.

When the tuning parameters for a complex model, including mechanistic considerations as well as empirical aspects, are adjusted for a very restrained database, they are only valid within this limited range. In this case little or no additional information regarding the underlying physic-chemical processes of the system can be gained, regardless of the level of detail of the model. Simple empirical correlations can often give equivalent results than a continuous lumping model within the restrained range of applicability. This work aims to show that a continuous lumping model can be validated for a large range of operating conditions (temperature, contact time, H2/HC), and feedstocks (residue content, nitrogen content).

Comprehensive reviews of the recent development in hydrocracking models have been compiled by Ancheyta and Elizade [START_REF] Ancheyta | Kinetic modeling of hydrocracking of heavy oil fractions: A review[END_REF][START_REF] Elizalde | Modeling the simultaneous hydrodesulfurization and hydrocracking of heavy residue oil by using the continuous kinetic lumping approach[END_REF]. A number of different approaches are used, depending on the available data and intended application of the model. Discrete lumping models, which use dedicated correlations for the reactivity of each lump, are among the simplest option. First developed by Qader and Hill [START_REF] Qader | Hydrocracking of Petroleum and Coal Oils[END_REF], these models have been widely used [START_REF] Ancheyta | Kinetic modeling of hydrocracking of heavy oil fractions: A review[END_REF][START_REF] Becker | Hydrotreatment modeling for a variety of VGO feedstocks: A continuous lumping approach[END_REF]. The lumps can represents specific cuts or pseudo components [START_REF] Krishna | Use of an axial-dispersion model for kinetic description of hydrocracking[END_REF]. Continuous lumping models offer a much finer resolution of the feed.

They consider the hydrocarbon mixture to be a continuous distribution along an internal coordinate (most commonly True Boiling Point, TBP) [START_REF] Ancheyta | Kinetic modeling of hydrocracking of heavy oil fractions: A review[END_REF][START_REF] Becker | Hydrotreatment modeling for a variety of VGO feedstocks: A continuous lumping approach[END_REF][START_REF] Elizalde | Modeling the simultaneous hydrodesulfurization and hydrocracking of heavy residue oil by using the continuous kinetic lumping approach[END_REF][START_REF] Lababidi | Modeling the hydrocracking kinetics of atmospheric residue in hydrotreating processes by the continuous lumping approach[END_REF][START_REF] Govindhakannan | On the construction of a continuous concentration-reactivity function for the continuum lumping approach[END_REF]. A more detailed description of the reaction network can be achieved by decomposition of the feed into a number of lumps, each distributed along TBP [START_REF] Becker | Hydrotreatment modeling for a variety of VGO feedstocks: A continuous lumping approach[END_REF]. Other modeling approaches include the use of neural networks [START_REF] Elkamel | MODELING THE HYDROCRACKING PROCESS USING ARTIFICIAL NEURAL NETWORKS[END_REF] and microkinetic modeling [START_REF] Narasimhan | Relumped singleevent microkinetic model for alkane hydrocracking on shape-selective catalysts: catalysis on ZSM-22 pore mouths, bridge acid sites and micropores[END_REF][START_REF] Martens | Kinetic Modeling of Paraffins Hydrocracking based upon Elementary Steps and the Single Event Concept[END_REF][START_REF] Guillaume | Single Event Kinetic Modelling without Explicit Generation of Large Networks: Application to Hydrocracking of Long Paraffins[END_REF][START_REF] Surla | Kinetic Modeling using the Single-Event Methodology: Application to the Isomerization of Light Paraffins[END_REF]. Using neural networks is an entirely empirical black-box approach, which is not based on an understanding of the underlying chemical kinetics. These powerful models can be used when large amounts of data are available. Microkinetic modeling, such as the single events approach, is based on the reconstruction of the feed into its individual molecules. A large reaction network with thousands of possible reaction pathways is then constructed. These models require a detailed understanding of the chemical kinetics and feed composition. A large amount of computational resources is also required.

The data fitting procedure and statistical analysis of model parameters for a continuous lumping model is presented in this work. The paper is structured as follows: The experimental setup of the pilot plant and the databases is presented in section 2. Details of the continuous lumping model used here are given in section 3. This is followed by a description of the parameter fitting procedure, as well as a model evaluation in sections 4 and 5. A statistical analysis of the model parameters is given in section 6.

Materials & Methods

Process Description

Industrial hydrocracking operations are generally a two-step process (R1 + R2). The second step (R2) is considered separately in this study.The main HCK reactor, with a zeolite catalyst follows an hydrotreatment reactor. The first reactor (R1) uses a catalyst designed to remove organic nitrogen and sulfur, hydrogenate aromatic compounds, and remove metallic impurities. These catalysts are much less susceptible to poisoning than the zeolite-based HCK catalysts. Nitrogen-and sulfur-containing compounds, as well as metals are known to inhibit zeolite catalysts [START_REF] Rana | A review of recent advances on process technologies for upgrading of heavy oils and residua[END_REF]. This allows the more fragile zeolite catalyst in the second reactor (R2) to perform the main hydrocracking reaction. The main purpose of the HCK catalyst (R2) is the breaking up of large hydrocarbon atoms contained in the VGO residue, in order to obtain the more valuable middle distillate cuts. Two types of tests were used in the calibration and evaluation databases: 1) pretreatment (R1) and hydrocracking (R2) in a single step, and; 2) separate pretreatment (R1). These two cases are illustrated in Figure 1 In the first case, the gases produced in R1 (H2S, NH3, hydrocarbon gasses) are carried over in the feed entering R2. The total mass entering R2, with respect to the feed of R1, is equal to the sum of the mass entering two-step process and the hydrogen consumption in R1.. The total yield of the two reactors is then equal to the mass of the feed (R1) and the hydrogen consumption in the two reactors (R1 and R2).

In the second case, the gases are separated from the effluent of the pretreatment reactor, which is analyzed and stored for subsequent use in a hydrocracking test. The feed entering R2 is therefore 100% liquid. Aniline and dimethyldisulfide (DMDS) additives are added in order to include the effect of NH3 and H2S on the zeolite catalyst. This mass is subsequently subtracted from the effluent, giving a total yield of sum of the mass of the feed (R2) and the hydrogen consumption (excluding H2 used up in the generation of H2S from DMDS and NH3 from Aniline). This setup allows the influence of feed nitrogen content to be studied independent of other feed characteristics.

The sulfur and nitrogen contents in the liquid stream entering R2 are subsequently referred to as S R2 and N R2 , respectively. The nitrogen and sulfur content entering R1 are referred to as N R1 and S R1 . In the case of a pre-treated feed, the equivalent nitrogen content due to aniline additive is used instead.

The Pilot Plant

The experimental runs presented in this study were performed in a pilot unit at IFP Energies Nouvelles, Solaize, France.

The hydrocracking step was performed on a commercial NiMoP-zeolite catalyst. Various commercial HDT catalysts were used for the pretreatment step. Total catalyst volume in both reactors is 50 cm 3 . The plant consists of a number of fixed beds, down-flow reactors, designed to mirror the operating conditions in industrial hydrocracking units. Unlike industrial units, which operate in adiabatic mode, the pilot plant operates in isothermal conditions. Temperature is controlled along the reactor. The units were run in continuous operation. Individual mass balances were taken for up to 12 hours, after temperatures, pressures, flow rates, and effluent properties were stabilized. A series of mass balances with different operating conditions were thus taken from each experimental run. Each mass balance corresponds to a single experimental point.

Analyses were performed on the feedstocks, the liquid and gaseous effluents, as well as a sample taken on the effluent of the pretreatment reactor (where applicable, see Figure 1 and Figure 2). The most relevant measurements to this study were the analysis of the feed sulfur and nitrogen contents, the partial pressures of NH3, H2S, and H2 gasses, as well as the simulated distillation.

The Calibration & Validation Databases

A calibration database, consisting of 52 mass balances, was compiled. Feed characteristics and operating condition were chosen to correspond to the ones encountered in industrial hydrocracking units. This database was used for the identification of the empirical parameters in the continuous lumping model. A validation database, consisting of an additional 22 mass balances was used to validate the model performance. The feed nitrogen and sulfur content of the two databases is shown in Figure 3. Different feeds of Iranian or Arabian origin were used at varying levels of pretreatment (R1). For about half (n=31) of the experimental runs, pretreatment and hydrocracking were performed in a single step, corresponding to Figure 1. The remaining runs (n=43) were performed in separate steps, corresponding to Figure 2.

The range of the main operating conditions, temperature (°C), hydrogen to hydrocarbon ratio (H2/HC, StdL/L), and Liquid Hourly Space Velocity (LHSV, h -1 ) is given in Figure 4. The residue (>370°C) cut in the feed of the zeolite catalytic bed (R2) and its conversion, X 370+ =100*(1-Y 370+,out /Y 370+,in ) are given in Figure 5 to give an overview of the performance of the experimental runs.

The range of the residual organic nitrogen in the feed is 0 -154 ppm. Four runs in the calibration database had contained no detectable organic nitrogen; these are not shown in Figure 3 because of the log scale. The sulfur content is not used in the modeling because it has been found to have negligible influence on hydrocracking reactions on a zeolite catalyst. The NH3 entering R2, due to the feed nitrogen content of the pretreatment reactor or the aniline additive was varied between 570 and 2800 ppm, with one run performed at a very high value of 5000 ppm. The nitrogen and sulfur content in the R1 and R2 feeds of the validation database are within the range of the validation database. The variations in the feeds must be taken into account by the model. Using a continuous rather than a discrete lumping approach allows the full distribution (SIMDIS) of the feed rather than a limited number of lumps to be used as the model input.

Most tests were performed at a standard ratio of H2 gas to hydrocarbon flowrates (H2/HC) of 1000 L/L (std) and a total pressure of 140 bar. Four runs in the calibration database were performed at 100 bar, while four runs in the validation database were performed at 160 bar. The temperature and LHSV were varied between 370 -401°C and 1.3 -4.2 h -1 in the calibration database. The validation tests generally fall within the same range, with the exception of four runs with a LHSV of 0.9 h -1 , and one test at a very low temperature of 362°C.

The residue (>370°C cut) of the feed of R2 falls in a relatively narrow range of between 67 -84 %, with two validation tests at around 50 %, corresponding to very high levels of pretreatment. The conversion (X 370+ ), on the other hand, covers the complete range that can be achieved in a hydrocracking reactor, i.e. up to 100 %.

Measurements for the composition of the residue (>370°C) cut of the calibration feeds (R2), in terms of paraffin (P), naphtene (N) and aromatic (A) compounds are available for the pretreated feeds (see Figure 2). This shows that the composition of the feeds used in this study does not vary very much. The residue cut of the feeds have an aromatics content of 6 137 7

Model Development

The Continuous Lumping Framework

The modeling framework considered the hydrocarbon mixture to be a continuous distribution of true boiling point (TBP). The reactivity of the components follows a curve described as a function of normalized TBP θ. Hydrocarbons below 8 carbon atoms (TBP ≈ 124°C) are assumed to undergo no further cracking, while cracking of very heavy hydrocarbons (>700°C) is considered to be constant and thus independent of TBP. The normalized TBP is defined as:

𝜃𝜃 = 𝑇𝑇𝑇𝑇𝑇𝑇 -𝑇𝑇 𝑚𝑚𝑚𝑚𝑚𝑚 𝑇𝑇 𝑚𝑚𝑚𝑚𝑚𝑚 -𝑇𝑇 𝑚𝑚𝑚𝑚𝑚𝑚 (1)
Where:

TBP True Boiling Point (°C) T min -7.6°C (TBP of C4H10) T max 700°C
It is important to note that the main cracking mechanism is β-scission of saturated hydrocarbon chains or cycles. βscissions on a straight chain result in two smaller molecules with lower boiling points. For endo-cyclic β-scissions of cycloalkanes (i.e. ring opening), a single, saturated alkane chain is formed. This results in a change in TBP, but no reduction in size of the molecule. However, a number of different reaction types occur during the hydrocracking process. Isomerization reactions change the branching degree of a hydrocarbon, which can increase or decrease the TBP of the reaction product. Aromatic compounds are hydrogenated, forming saturated cycloalkanes (naphtenes). Hydrogenation of aromatic ringsis also accompanied by a decrease in boiling point.

These reactions are not explicitly modeled in the continuous lumping model presented here. The set of reactions, causing a decrease in the TBP of the mixture are lumped together to give an apparent reaction rate. A more detailed, multifamily approach, which treats the reaction pathways between aromatic, naphtene, and paraffins explicitly, has been presented by Becker et al. [START_REF] Becker | Hydrotreatment modeling for a variety of VGO feedstocks: A continuous lumping approach[END_REF]. Fitting the parameters of this model requires detailed knowledge of the feed and effluent composition.. Such analyses are rarely available in an industrial context. A model capable of predicting the composition of the effluent is therefore of limited usefulness for the hydrocracking units design and optimization. However, they are useful in a research and development context. They allow the characterization of the chemical reaction pathways and activity of a catalyst. This can lead to a better understanding of the underlying phenomena which can consequently be used to guide future developments.

The model presented here has been developed to capture the effect of the main operating parameters and feed conditions in order to predict the yield structure of the hydrocracker.

The Model Equations

The feed simulated distillation (SIMDIS is a cumulative mass distribution. The associate continuous distribution, c(θ,t), can be split up into a number of intervals i, with Δθ i = [𝜃𝜃 𝑚𝑚-, 𝜃𝜃 𝑚𝑚+ ], with 𝜃𝜃 𝑚𝑚-< 𝜃𝜃 𝑚𝑚+ . The discrete concentration in each interval, C i , is given by: .

𝐶𝐶 𝑚𝑚 = � 𝑐𝑐(𝜃𝜃, 𝑡𝑡)𝑑𝑑𝜃𝜃 𝜃𝜃 𝑖𝑖-1 𝜃𝜃 𝑖𝑖-1 (2)
Each interval is represented by θ i . The mass balance for each interval i is:

Δ𝐶𝐶 𝑚𝑚 Δ𝑡𝑡 = -𝑅𝑅 𝑚𝑚,𝑜𝑜𝑜𝑜𝑜𝑜 + 𝑅𝑅 𝑚𝑚,𝑚𝑚𝑚𝑚 (3)
The mass being eliminated from the interval i, R i,iout is due to molecules being cracked. This results in a lower TBP. The mass R i,in being added to the interval i, is due to the cracking of heavier hydrocarbons (i.e. in classes with θ j > θ i ). The mass of hydrocarbons resulting from a cracking reaction, which fall in the Δθ i interval, i.e. the R i,out term, is added to the mass

C i .
A set of partial differential equations for the hydrocracking reaction can be constructed for the continuous mass distribution c(θ,t), given in equation [START_REF] Motaghi | Slurry-phase hydrocracking-possible solution to refining margins: Opportunity crudes require more hydrogen addition to upgrade orphan product streams into higher-value "clean" products[END_REF] . Here, the continuous elimination and creation terms are denoted r out and r in respectively. The yield distribution, g(θ,θ * ), and reactivity function, f X (θ), do not change along the reactor (for isothermal operation)

𝜕𝜕𝑐𝑐(𝜃𝜃

and are therefore calculated a priori.

A gamma function was chosen for the reactivity function and a beta distribution was chosen for the yield distribution.

These two functions are very flexible regarding the shapes they can take [START_REF] Becker | Hydrotreatment modeling for a variety of VGO feedstocks: A continuous lumping approach[END_REF]. They were chosen over the mode common Gaussian and exponential distribution to give additional degrees of freedom to the model. The functions are given in equations ( 4) and ( 5) respectively Both functions take two tuning parameters. The yield distribution is defined in terms of carbon numbers (n), rather than TBP (θ) in order to ensure mass conservation without an additional normalization term. The model equation (3), however, is defined in terms of TBP (θ). The transformation of carbon number to TBP, 𝑛𝑛 = 𝑓𝑓(𝜃𝜃), is performed according to the TBP of normal paraffins.

𝑓𝑓(𝜃𝜃) = 𝜃𝜃 𝐴𝐴 𝑒𝑒 -𝐵𝐵𝜃𝜃 (6)
An Arrhenius term was included in the pre-factor (k T ) to account for the effect of temperature variations.

𝑘𝑘 𝑇𝑇 = 𝑘𝑘 0 exp �- 𝐸𝐸 𝑚𝑚 𝑅𝑅 𝑔𝑔 � 1 𝑇𝑇 - 1 𝑇𝑇 𝑟𝑟𝑟𝑟𝑟𝑟 �� (8)
With:

k 0 Pre-factor,model parameter

E a

Activation energy,model parameter (in kJ mol -1 ) T Reaction temperature (in K)

T ref Reference temperature (in K) 11 R g Gas constant (8.314 J K -1 mol -1 )
The continuous expression (equation 3) is discretized at equally spaced TBP intervals of 5°C. The resulting set of ODEs is solved by the well-known LSODE solver. The model was implemented in object-oriented FORTRAN.

Hydrogen Consumption

The total mass flowrate of the effluent is higher than the mass flowrate of the (HCK) feed because of the hydrogen added in process. Estimation of the hydrogen consumption is necessary in order to determine the total yield of the hydrocracking reaction. Two hydrogen atoms (i.e. one H2) are added in each β-scission reaction. The total number of cracking reactions is estimated via a correlation of the feed and effluent densities, i.e. 𝐻𝐻 2,𝛽𝛽 = 𝑓𝑓�𝑑𝑑 154,𝑟𝑟𝑟𝑟𝑟𝑟𝑓𝑓 , 𝑑𝑑 154,𝑟𝑟𝑟𝑟𝑟𝑟 �. This correlation is based on the n-d-m method (ASTM D-3238 [START_REF] Astm D | Standard Test Method for Calculation of Carbon Distribution and Structural Group Analysis of Petroleum Oils by the n d M Method n[END_REF]) is used.

Furthermore, three H2 molecules are added for each full hydrogenation of an aromatic ring: 𝐻𝐻 2,𝐴𝐴 = 3 * (𝐶𝐶 𝐴𝐴,𝑟𝑟𝑟𝑟𝑟𝑟𝑓𝑓 -𝐶𝐶 𝐴𝐴,𝑟𝑟𝑟𝑟𝑟𝑟 ).

Measurements of the total amount of aromatic carbon atoms in the feed are available, and the extent of the hydrodearomatisation reaction is estimated by a dedicated model not detailed here. The hydrogen added by the hydrodenitrogenation and hydrodesulfurization reactions is also taken into account.

Inhibition Terms

Zeolite catalysts are known to be inhibited by the presence of NH3 gas [START_REF]Catalyse acido-basique[END_REF]. An inhibition term (I N ), based on the nitrogen content in the feed of R1 (N R1 ) was therefore introduced in the model. The nitrogen in the feed of R1 is mostly transformed to NH3, which is carried over to R2. The inhibition term is therefore derived from an estimation of the partial pressure of NH3 gas in the reactor, given in equation [START_REF] Martens | Attempts to Rationalize the Distribution of Hydrocracking Products. I Qualitative Description of the Primary Hydrocracking Modes of Long Chain Paraffins in open[END_REF].

1 ppNH 3 ∝ � 𝐻𝐻 2 𝐻𝐻 𝐶𝐶 𝑇𝑇 𝑜𝑜𝑜𝑜𝑜𝑜 𝑁𝑁 𝑅𝑅1 � (9)
The residual nitrogen content in the feed of R2 (N R2 ) does also exhibit an inhibition effect on the zeolite catalyst. The nitrogen containing organic compounds adsorb on the acidic sites of the catalyst. These sites are therefore no longer available to perform hydrocracking reactions, until the nitrogen has been converted into NH3 gas. The nitrogen containing compounds remaining after a hydrotreatment step are generally very difficult to break down. The complete inhibition term used in this model is given in equation [START_REF] Henry | Vacuum gas oil hydrocracking performance of bifunctional Mo/Y zeolite catalysts in a semi-batch reactor[END_REF]. A total of four tuning parameters allow the adjustment of this term to experimental data.

�1 + 𝛼𝛼 𝑁𝑁𝑅𝑅2 𝑁𝑁 𝑅𝑅2 𝛽𝛽 𝑁𝑁𝑁𝑁2 � � ⎣ ⎢ ⎢ ⎡ 1 1 + 𝛼𝛼 𝑁𝑁𝑅𝑅1 � 𝐻𝐻 2 𝐻𝐻 𝐶𝐶 𝑇𝑇 𝑜𝑜𝑜𝑜𝑜𝑜 𝑁𝑁 𝑅𝑅1 � 𝛽𝛽 𝑁𝑁𝑁𝑁1 ⎦ ⎥ ⎥ ⎤ (10)
With:

H2H C

Volumetric ratio of H2 to hydrocarbon feed (in std L/L )

P tot

Total pressure in R2 (in bar)

N R1
Organic nitrogen in feed of R1 (in ppm)

α NR2 , β NR2 Model parameters N R2
Organic nitrogen in feed of R2 (in ppm) α NR2 , β NR2 Model parameters

Generation of Hydrocarbon Gas

The generation of hydrocarbons smaller than C5 is treated separately in this model. It has been observed that the process of generating these light gases is rather different and cannot be modeled by the same continuous yield distribution used for the simulation of the cracking of heavier hydrocarbons. A separate term, g(θ,GAS), was therefore introduced in the yield distribution.

𝑔𝑔(𝜃𝜃, 𝐺𝐺𝐺𝐺𝐺𝐺) = 𝛼𝛼 𝑔𝑔𝑚𝑚𝑔𝑔 𝜃𝜃 𝛽𝛽 𝑔𝑔𝑔𝑔𝑔𝑔 𝑒𝑒 -𝜃𝜃 (11)
This term, shown in equation [START_REF] Burnens | Hydrocracking reaction pathways of 2,6,10,14tetramethylpentadecane model molecule on bifunctional silica-alumina and ultrastable Y zeolite catalysts[END_REF], uses a separate gamma function to estimate the total yield of C1-C4 gas by cracking of a hydrocarbon with TBP θ. The additional model parameters (α gas and β gas ) give further flexibility to the model.

The redistribution between C3 and C4 gasses are calculated using a dedicated correlation, which has been fitted separately using a separate calibration database. Zeolite catalysts, such as the one used in this study, generate negligible amounts of C1 and C2 gas. This is because the β-scission reaction mechanism cannot crack a hydrocarbon chain at the ends.

Parameter Identification

The continuous lumping model presented in section 2.1 contains 12 tuning parameters. These parameters were identified using the calibration database presented in section 3.3. An objective function for a set of parameters (β) was constructed, defined in equation ( 9), from the square difference between the experimental measurements (y) and the simulation results 𝑓𝑓(𝑥𝑥, 𝛽𝛽). Weights (w) were applied different components of the objective function. The yield structure (standard cuts & gasses C1-C4), the SIMDIS, and the total yield (100% + conso H2), as well as the SIMDIS (9 points), were used in the objective function, according to following weights (w):

• Naphtha (< 150°C cut) 3 • Kerosene (150-250°C cut) 6 • Gasoil (250-370°C cut) 6 • Middle distillates (150-370°C cut) 10 • Residue (> 370°C cut) 10 • SIMDIS (10 -90%) 0.01 • Gas yield (C1-C4) 10 • Total yield 1
It is important to note that the error in SIMDIS is expressed in terms of °C, while the other errors are expressed in terms of percentages. The former are an order of magnitude higher than the latter for a comparable fit of the model. A low weight was therefore chosen in order to obtain the same order of magnitude for all observables.

The DN2FB FORTRAN optimization routine, part of the PORT library (Bell Labs) was used for parameter identification [START_REF] Dennis | Algorithm 573: NL2SOL---An Adaptive Nonlinear Least-Squares Algorithm [E4[END_REF]. This routine uses a Levenberg-Marquardt least squares algorithm to find the gradients of the objective function. The target function is highly non-linear. The optimized set of parameters was found to be highly dependent on the initial parameters, which implies convergence to local minima, rather than global ones in some cases. Optimizations were therefore performed for a set of 96 randomly chosen initial parameters. Optimizations were performed on the highperformance cluster at IFP Energies Nouvelles, Solaize, France. The parameters of the case which converged to the lowest value of the objective function was retained.

This procedure does not guarantee convergence to the global minimum, it does, however, at the very least provide a good local minimum. The use of more advanced, global algorithms (i.e. genetic / global response surface algorithms) is far more complex and time consuming than the gradient-based algorithm used in this study. This might be justified in a pure academic research context.The approach used in this study is better suited for practical applications, where a the model parameters for a new catalyst needs to be identified rapidly. The comparison of different optimization algorithms is performed at IPFEN and is the object of separate bodies of work.

Model Validation

The objective of the construction of a continuous lumping model for a hydrocracking process is the accurate estimation of yield structure and simulated distillation for a range of operating conditions. Estimation of the yield structure, in terms of standard cuts, is necessary in order to perform process design and optimization. It is common to design a hydrocracking process in order to maximize the conversion of the residue cut (X 370+ ) and the yield of the middle distillates (150 -250°C cut). The production of hydrocarbon gases (C1-C4) is generally undesirable [START_REF] Rana | A review of recent advances on process technologies for upgrading of heavy oils and residua[END_REF].

The temperature of an industrial hydrocracking reactor is often adjusted to give a desired conversion (X 370+ ). The relationship between X 370+ and yield structure is therefore of particular interest. For a given feed and catalyst, this relationship is relatively independent of reactor temperature. Simulations are run in two different modes:

1) iso-T: the real (measured) reactor temperature is used;

2) iso-RS: the temperature is adjusted such that the simulated residue content equals the experimental measured value.

Good estimation of the SIMDIS of the effluent is necessary in order to calculate a number of product qualities for the individual cuts (e.g. d154, octane number, viscosity index, cetane index). This is done by dedicated correlations, based on the modeled SIMDIS.

Parity graphs, showing the agreement of simulation results and experimental measurements have been constructed for iso-T and iso-RS modes.

Iso-T mode

The results for residue conversion (X 370+ ) and middle distillate yield in iso-T mode are shown in Figures 8 and9.

The residue conversion is well predicted for the entire range, i.e. 15 to 100%. Some data points, particularly for the validation database and/or for very high X 370+ , tend to be somewhat under-predicted. The middle distillate yield (150-370°C) is very well predicted in iso-T mode. All points, with the exception of one point in the validation database, fall within an error band of ±5 %. A summary of the mean absolute errors for a number of indicators is shown in Table 1.

The mean absolute errors for X 370+ and yield structure for both, the calibration and validation database are very good for experimental runs with X 370+ < 90%. The errors are slightly higher but remain acceptable for the validation points. The points with very high residue conversion are less well predicted by the model. It is, however, important to note that industrial hydrocracking units rarely operate at conversions higher than 90%. These points are therefore of somewhat limited interest.

The possible reasons for this behavior are that the calibration database contained only two data points with X 370+ > 90%.

The shape of the reactivity function in this region is therefore not sufficiently sensitized. The other reason is related to the fact that the SIMDIS becomes very steep around 370°C. This amplifies the error of residue ( > 370°C cut) yield, i.e. the point at which the SIMDIS and intersects the 370°C line. The mean absolute errors of hydrogen consumption and total yield show that the model ensures mass conservation and predicts the correct mass flowrates.

Iso-RS mode

A more detailed analysis is provided for the simulation results in iso-T mode. The simulation results for residue (>370°C cut) yield and thus X 370+ , are fitted to experimental measurements in this mode. The yield structure of the model is expected to improve with respect to iso-T mode. Figure 10 shows the middle distillate (150-370°C cut) yield. Kerosene (150-250°C cut) and gasoil (250-370°C cut), which make up the middle distillate are shown in Figures 11 and12 respectively.

The middle distillate (150-370°C cut) yield is better predicted than in iso-T mode for points with X 370+ < 90%. The majority of points fall within the ±2 % error band. The high conversion points, on the other hand, are less well predicted in iso-RS mode than in iso-T mode. The parity graphs for the gasoil and kerosene yields show that the error in middle distillate comes mainly from the 250-370°C cut. Kerosene (150-250°C cut) yield is well predicted, even for elevated residue conversion. The gasoil (250-370°C cut) yield is poorly predicted for experimental runs with X 370+ >90%, even though most points fall within a very narrow range (~20-25%).

The parity graph for total hydrocarbon gas (C1-C4) shows that the continuous lumping model provides very good simulation results, even for high conversion (X 370+ > 90%). Only three points of the validation database fall outside the error band of ±1 %. A summary of the mean absolute errors for a number of indicators is shown in Table 2.

The ΔT values reported in Table 2 correspond to the mean temperature difference that had to be applied in order to obtain the desired X 370+ . As in iso-T mode, the results for the validation database are slightly degraded with respect to the results obtained for the calibration points. A ΔT of 2.0°C and an error of the yield structure of slightly more than 1 % are well within the accuracy that can be expected from a hydrocracking model.

influence on the model behavior. The parameters of the gas distribution (equation 9) have the least influence on the objective function.

The correlation between the individual parameters can be determined from the symmetry in the Eigenvectors. This is shown in Table 4, higher values strong correlations. This table shows that the pre-exponential factor (k 0 ) is strongly correlated with the inhibition term for NR1. The two parameters in the NR2 inhibition term (α NR2 and β NR2 ) are strongly correlated. The same can be said regarding the parameters governing the yield distribution for the gasses (Equation 11).

Conclusions

A continuous lumping model was successfully applied to hydrocracking of VGO feedstock on zeolite catalysts. The second reactor in a two-stage hydrocracking process was simulated. The model contains inhibition terms for NH3 carried over from the pretreatment reactor (R1), as well as the residual organic nitrogen in the feed of the hydrocracking reactor (R2). A gamma function was used for the continuous reactivity function and a beta distribution was used for the yield distribution.

The 12 model parameters were identified using a calibration database with 52 data points. Parameter identification was performed using 96 sets of randomized initial parameters. The model was then used to simulate a validation database with 22 additional data points. Good simulation results, in terms of X 370+ and yield structure were obtained. The model performed less reliably in the case of high conversion (X 370+ > 90%).

A statistical parameters analysis, based on the Eigenvectors of the Hessian around the optimum was presented. This allowed the parameters with the largest impact on the objective function to be identified. An analysis of the correlation between the parameters was also presented. The eigenvector analysis was used to show that the nitrogen inhibition terms have (NR1 and NR2) have a dominating effect on the target function. This makes it clear that the inclusion of these terms is absolutely necessary to correctly simulate hydrocracking of feeds with different nitrogen content and different levels of pre-treatment. Furthermore, a strong correlation between the nitrogen (NR1) inhibit terms and the pre-exponential factor has been shown.

4 )

 4 , 𝑡𝑡) 𝜕𝜕𝑡𝑡 = -𝑟𝑟 𝑜𝑜𝑜𝑜𝑜𝑜 (𝜃𝜃, 𝑡𝑡) + 𝑟𝑟 𝑚𝑚𝑚𝑚 (𝜃𝜃, 𝑡𝑡) (An expansion of the elimination (r out ) and generation (r in ) terms leads to: 𝜕𝜕𝑐𝑐(𝜃𝜃, 𝑡𝑡) 𝜕𝜕𝑡𝑡 = -𝑘𝑘 T 𝑓𝑓 𝑋𝑋 (𝜃𝜃)𝐼𝐼 𝑁𝑁 (𝑡𝑡)𝐶𝐶(𝜃𝜃, 𝑡𝑡) + � 𝑘𝑘 T 𝑓𝑓 𝑋𝑋 (𝜃𝜃)𝐼𝐼 𝑁𝑁 (𝑡𝑡)𝐶𝐶(𝜃𝜃, 𝑡𝑡)𝑔𝑔(𝜃𝜃, 𝜃𝜃 * )

n

  Number of carbon atoms in compound undergoing cracking n * Number of carbon atoms in cracking product α y , β y Model parameters for yield distribution A, B Model parameters for reactivity

Statistical Analysis of Model Parameters

Statistics were calculated by an in-house code to quantify the relative importance and correlations between the parameters. Standard non-linear optimization routines, such as DN2FB, attempt to determine the optimal set of parameters for the best local minimum of the least-squares objective function. This implies a zero gradient hat around the optimal set of parameters β 0 :

The Taylor expansion of the objective function (J), around the optimal parameters (β 0 ), given the Hessian (H) gives:

The Eigenvector (v), with the associated Eigenvalue (μ) of the Hessian (H) is defined as:

A perturbation of β in the direction of the Eigenvector (v) gives then, by definition:

𝐽𝐽(𝑥𝑥, 𝛽𝛽) = 𝐽𝐽(𝑥𝑥, 𝛽𝛽 0 ) + 1 2 𝛼𝛼𝐻𝐻′𝐻𝐻𝛼𝛼𝐻𝐻

( 18)

𝐽𝐽(𝑥𝑥, 𝛽𝛽) = 𝐽𝐽(𝑥𝑥, 𝛽𝛽 0 ) + 1 2 𝜇𝜇𝛼𝛼 2 ||𝐻𝐻||²

( 19)

An analysis of the components of the Eigenvector and the magnitude of the associated eigenvalues allows the relative impact of variations of the individual parameters on the target function to be quantified. A small eigenvalue implies a small impact of a variation of the parameters in the direction of the associated Eigenvector; inversely, a large eigenvalue implies a large impact of a perturbation of the parameters in the direction of the associated eigenvector.

The Eigenvectors with their associated Eigenvalues are shown in Table 1. The dominant parameter in the first five Eigenvectors are highlighted to show the parameter with the largest impact on the objective function. This shows that the model is very sensitive to the exponents of the nitrogen inhibition term (β NR1 and β NR2 in equation 8). The parameters of the reactivity function (equation 4), the activation energy (E a ), and the yield distribution (Equation 7) also have a large