Memoria Para

Optar Al

DOCTOR EN CIENCIAS DE LA INGENIERÍA MENCIÓN MODELACIÓN Grado De

Alejandro Maximiliano

Rojas Palma

Prof Guía: Héctor

Ramírez Cabrera

Alain Rapaport

Estudio De

Algunos Problemas

E N Modelación

Y Optimización

D E Bioprocesos

Résumé De

POUR OBTENIR L A Thèse

L E Degre

D E Docteur

Alejandro Maximiliano Rojas

Palma Date

Ramírez Héctor

Alain Cabrera

Etude Rapaport

Problèmes De Certains

La Dans

L ' Modélisation

Des Optimisation

Bioprocédés

Rapaport

RESUMEN DE LA

Keywords: Bioprocesos, optimización, modelación matemática, sistemas compartimentados, crecimiento de microalgas, controlabilidad Bioprocédés, optimisation, modélisation mathématique, systèmes compartimentaux, croissance des microalgues, contrôlabilité Bioprocess, optimization, mathematical modelling, compartmental system, microalgae growth, controllability

El propósito de esta tesis es

Agradecimientos

En primer lugar quiero agradecer a mi director de tesis, Héctor Ramírez Cabrera, por su constante apoyo a largo de esta tesis, así como también agradecer el haber contado con el apoyo del profesor Alain Rapaport, sobre todo por sus indicaciones durante mis pasantías de investigación en el Campus INRA-Supagro en Montpellier, Francia.

Agradezco también a los profesores David Jeison (UFRO) y Francis Mairet (INRIA-Sophia Antipolis) por sus sugerencias y comentarios, así como por las participaciones en reuniones de trabajo y talleres relacionados al grupo BIONATURE al cual pertenecí en calidad de alumno tesista de doctorado. A Gonzalo Ruiz † , Andrés Donoso, Fabio Carrera y al grupo de investigación de la escuela de Ingeniería Bioquímica de la PUCV por sus consejos, en particular en lo que respecta a la interpretación de resultados y explicaciones bioquímicas de los procesos estudiados. También debo mencionar a Matthieu Sebbah (UTFSM), pues sin su ayuda no hubiese podido comunicarme en mi primera pasantía, a Mario Veruete (UMontpellier) por su ayuda en todo lo referido a la gestión en Francia durante todas mis estadias y a Victor Riquelme (UChile) por sus aportes en aspectos matemáticos en los que discutimos con frecuencia, llegando a buenos resultados.

Agradezco a CONICYT por financiar mis estudios a través de la beca de doctorado nacional y de la beca de tesis de doctorado en la empresa. Al Departamento de Postgrado y Postítulo de la Vicerrectoría de Asuntos Académicos de la Universidad de Chile y al Instituto Francés de Chile de la Embajada de Francia por financiar algunas de mis pasantías a traves del programa de becas de movilidad doctoral de estadías cortas de investigación en Francia para estudiantes tesistas. También agradezco a INRIA Chile por permitirme ser parte del proyecto BIONATURE y con ello participar en diversos workshops organizados en el país, además de ser patrocinantes en mi proyecto de tesis en la empresa. Su ayuda pecuniaria fue fundamental.

Finalmente quisiera agradecer a mi familia por todo su apoyo durante estos años, en particular a mi esposa Alejandra, a mis padres Nino y María por ser vivos ejemplos de entereza, perseverancia y modelos a seguir. Finalmente a mi pequeña Consuelito, una luz en mi vida, a quien le dedico este escrito. [START_REF] Ballyk | Global analysis of competition for perfectly substitutable resources with linear response[END_REF] The standard n-vessel gradostat. The left vessel labeled R is a reservoir containing nutrient at concentration S 0 , C is an overflow vessel, and D denotes the dilution rate. All vessels have the same volume (source [START_REF] Smith | The Theory of the Chemostat[END_REF]). . . . . . . . . . . . . . . . . . . . 1.9 Irreducible "dead-end" gradostat. Note that the inflow to each vessel balances the outflow (source [START_REF] Smith | The Theory of the Chemostat[END_REF]). . . [START_REF] Ballyk | Exploitative Competition in the Chemostat for Two Perfectly Substitutable Resources[END_REF] The figure shows the results of the table 2.4. It may be noted that the higher the incidencial light, higher the final biomass concentration, while the terminal batch time decreases and the mean volumetric productivity increases, the initial concentration needed to achieve optimal productivity decreases slightly. . . . . . . . . . . . . . . . 2.10 Linear regression and zero intercept for the estimation of the light attenuation parameter a. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 Approximation of functions µ 1 (.) and µ 2 (.) by nonlinear least squares method. . . . .

3.1

In this setting μn = 2, μs = 1, K n = 1, K s = 0.12, K x = 0.5, s in = 1, D = 0.5, k 1 = 1 and considering the growth functions in the application section, the graph show the existence of the three equilibrium points that was proven in lemma above. . . . . 3.2 Graph of the function µ 3 (.) for k 2 = 1 and the parameter values in table 3 3.5 (from [1]) with a pond depth h = 0.1 m (left) and h = 0.5 m (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Solutions of the system (3.6) with initial condition (x 1 (0), x 2 (0), s 1 (0), s 2 (0)) = [START_REF] Berman | Nonnegative Matrices in the Mathematical Sciences[END_REF][START_REF] Gelhar | A critical review of data on field-scale dispersion in aquifers[END_REF][START_REF] Dutta | Fundamentals of biochemical engineering[END_REF][START_REF] Hsu | A survey of constructing Lyapunov functions for mathematical models in population biology[END_REF] and a pond depth h = 0.1 m, considering a dilution rate D = 0.1 (left) and D = 0.4 (right). In the left figure the equilibrium E a = (0, 1496, 0.045, 0) is globally stable, however in the right figure the positive equilibrium E c = (11.2461, 1034, 0.388, 36.88) is globally stable. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Solutions of the system (3.6) with initial condition (x 1 (0), x 2 (0), s 1 (0), s 2 (0)) = [START_REF] Berman | Nonnegative Matrices in the Mathematical Sciences[END_REF][START_REF] Gelhar | A critical review of data on field-scale dispersion in aquifers[END_REF][START_REF] Dutta | Fundamentals of biochemical engineering[END_REF][START_REF] Hsu | A survey of constructing Lyapunov functions for mathematical models in population biology[END_REF] and a pond depth h = 0.5 m, considering a dilution rate D = 0.9 (left) and D = 1.

General Introduction 1.1 Thesis description and motivations

At present there is a tendency to develop efficient industrial processes from economic and environmental point of view. In this sense, bioprocesses, which are defined as fermentation processes involving the decomposition of glucose contained in substrates by enzymes produced by microorganisms, are more interesting due to the production of several products (antibiotics, foods, beverages, enzymes, biofuels , etc.) in many productive areas of industry.

However, the mathematical modelling of these processes is valuable insofar as it answers theoretical questions concerning the kinetics of the involved reactions and optimizes the process itself under certain predetermined criteria, such as minimizing production time or overall costs or maximizing the biomass resulting from the bioprocess. This is the motivation for the study of the problems raised in this manuscript. Naturally, it is not possible to cover the broad spectrum of problems in this area; however, the purpose is to show a methodology of analysis that responds to the objectives and leads to discussions about the solutions obtained.

The focus of this manuscript is the proposal and analysis of three problems that arise from bioprocess engineering. The first problem involves optimal microalgae productivity in bioreactors with light incidence. Light has a positive influence on microalgae growth; thus, productivity is better in cultures with light incidence. A management objective in the industrial production of microalgae is productivity maximization over time. In this sense, in [START_REF] Grognard | Optimization of a photobioreactor biomass production using natural light[END_REF], the authors address the optimization of biomass long-term productivity within the framework of microalgal biomass production in photobioreactors in continuous operation mode under the influence of day/night cycles. They propose a simple bioreactor model accounting for light attenuation in the reactor due to biomass density and obtain a control law that optimizes productivity over a single day through the application of Pontryagin's maximum principle, with the dilution rate being the control. They conclude that because of the day-night constraint, the productivity rate cannot be as high as it could have been without it. However, when the maximal growth rate is sufficiently larger than the respiration rate, we manage to have a temporary phase where the productivity rate is at or near this level. The maximal harvesting at the end of the light phase and at the beginning of the dark phase minimizes the biomass during the dark phase and, consequently, the net respiration. If the maximal growth rate is very large, the optimal solution consists in constantly applying maximal control because the biomass that is built-up in the light phase needs to be harvested even during the night.

Naturally, a question arises: Is it possible to make a comparison regarding the behavior and optimization results considering a batch operation mode? The first substantial difference between the processes is related to the dilution rate. In continuous operation mode, as a chemostat model [START_REF] Smith | The Theory of the Chemostat[END_REF], the dilution rate is a parameter that can be considered a control variable. However, in a batch process, dilution occurs in the final time of the process because batch mode is a discontinuous process. Thus, from a mathematical point of view, an optimal control approach is not possible in this case. In the batch process there are two decision variables: the initial condition of biomass concentration of each species and the final time of the operation. The idea is to formulate a nonlinear optimization problem with the objective of maximizing the mean volumetric productivity. Following the same assumptions given in [START_REF] Grognard | Optimization of a photobioreactor biomass production using natural light[END_REF], the model of two species, substrate concentration and microalgae biomass concentration, in a batch culture is reduced to a first-order nonlinear ordinary differential equation with a discontinuous right-hand side [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides: Control Systems[END_REF] that represents the biomass concentration growth of microalgae under light incidence. The discontinuity of this equation is caused by the influence of day/night cycles that are given in the model as step functions (see appendices A and B in chapter 2). The solution of this equation is a continuous function but not differentiable in all the domains; where the solution is part of the objective function. This implies that the optimization problem formulated is non smooth and the classical gradient approach is not possible to apply. To solve this problem or at least give the necessary optimality conditions, some results of nonsmooth optimization theory were applied. The analysis, main results and discussion of this problem are given in chapter 2.

In the second problem of this thesis, we try to simplify a more complex model of microalgae growth in a chemostat in competition with oxidizing nitrite bacteria and oxidizing ammonia bacteria and various substrates namely, total inorganic carbon, ammonium, nitrite, nitrate and nitrous acid where nitrification is represented by a two-step bioreaction. The growth functions considered for each species take into account different inhibition and limitation terms. The mass balance equations form a high-dimensional system of nonlinear ordinary differential equations. This fact causes the mathematical analysis of the system to be overcomplicated and even unsuccessful. Although it is possible to apply a numerical approach to obtain the relevant information of dynamical behavior from a quantitative point of view, we expect to have an idea of the global behavior based on qualitative analysis.

The simplifying model is a four-dimensional nonlinear system that represents the dynamics of two species, microalgae and nitrifying bacteria, in competition for nitrogen present as ammonium and nitrate produced by nitrification in a continuous process (chemostat). The nitrification is represented by a one-step bioreaction in this case. The model that we formulate and the stability results, present several features that we believe to potentially be of interest. In this chemostat model, an intra-specific competition phenomenon is considered based on density dependent growth functions and also cross-feeding; i.e., the nitrate produced by the nitrifiers can be con-sumed by the algae. We show that the coexistence of species is possible in the chemostat, although their stability may vary depending on some parameter values. It is feasible to reduce the system using the theory of asymptotically autonomous systems, and the limiting system obtained can be considered as a perturbation of a system of two species competing for a substrate. Thus, we use a result of non-vanishing perturbed systems to obtain a strong stability theorem for equilibrium coexistence. Stability analysis, discussion and an application using realistic parameter values are provided in chapter 3.

The third problem studied in this manuscript is related to linear compartmental systems. In geoscience, models of fractured porous media are often described as a mobile zone driven by advection, and one or several immobile zones are directly or indirectly connected to the mobile zone by diffusion terms. We believe that these models are also relevant to describing flows in soil or in porous media such as biofilms. To model a flow process, it is possible to define different configurations between these zones. Example of these are the classical connections in series and parallel, both in non-interacting systems [START_REF] Coulson | Chemical and Biochemical Reactors and Process Control[END_REF]. The models MINC (Multiple INteracting Continua) and MRMT (Multiple Rate Mass Transfer) are extensively used in transport phenomena. For instance, MINC is used to model flow in fractured media [START_REF] Warren | The behavior of naturally fractured reservoirs[END_REF], and MRMT is used to model the transport of solutes both passive and reactive [START_REF] Haggerty | Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity[END_REF].

Recently, these models were studied in [START_REF] De Dreuzy | Influence of porosity structures on mixing-induced reactivity at chemical equilibrium in mobile/immobile multirate mass transfer MRMT and multiple interacting continua MINC models[END_REF]. Specifically, the authors are interested in the relation between these configurations. In a previous work (see appendix A in chapter 4) a transfer function approach is used to show that these two schemes are often used in the literature: MINC-where the diffusive compartments are connected in series-and MRMT-where the diffusive compartments are connected in a star around the mobile zone-are equivalent input-output representations and provide formulas (up to three compartments) to pass from one representation to another. This result means that one can simply choose the most convenient approach when dealing with control or optimization without loss of generality. In this thesis, we prove this equivalence in the general case of n compartments. For this, we analyze the equivalence between MINC/MRMT structures and a general class of input-output linear systems that represent solute transport and diffusion between different compartments. We indeed show that any such systems can be equivalently represented by any of the particular two matrix structures: MRMT or MINC (but using different volumes and diffusive transfer terms). Moreover, we propose explicit methods to build equivalent MRMT or MINC matrices from any given structure, on the condition that it fulfills a controllability assumption. As shown in the examples, this controllability condition is not provided by the single irreducibility assumption on the network structure. The results and examples mentioned above can be found in chapter 4.

Some important notions in bioprocess engineering

Bioprocess engineering is the application of engineering principles to design, develop, and analyze processes using biocatalysts. These processes may result in the formation of desirable compounds or in the destruction of hazardous substances. The tools of the engineer, particularly the chemical engineer, are essential to the successful exploitation of bioprocesses [START_REF] Shuler | Bioprocess Engineering: Basic Concepts[END_REF].

Biotechnology usually implies the use or development of methods of direct genetic manipulation for a socially desirable goal. Such goals might include the production of a particular chemical, but they may also involve the production of better plants or seeds or gene therapy or the use of specially designed organisms to degrade waste. The key element for many users is the use of sophisticated techniques outside the cell for genetic manipulation. Others interpret biotechnology in a much broader sense and equate it with applied biology; they may include engineering as a subcomponent of biotechnology [START_REF] Fitch | An Engineering Introduction to Biotechnology[END_REF].

Many words have been used to describe engineers working with biotechnology. Bioengineering is a broad title and would include work on medical and agricultural systems; its practitioners include agricultural, electrical, mechanical, industrial, environmental and chemical engineers, as well as others. Biological engineering is similar but emphasizes applications to plants and animals. Biochemical engineering has usually meant the extension of chemical engineering principles to systems using a biological catalyst to bring about desired chemical transformations. It is often subdivided into bioreaction engineering and bioseparations. Biomedical engineering has been considered to be totally separate from biochemical engineering, although the boundary between the two is increasingly vague, particularly in the areas of cell surface receptors and animal cell culture.

There is a difference between bioprocess engineering and biochemical engineering. In addition to chemical engineering, bioprocess engineering would include the work of mechanical, electrical, and industrial engineers who apply the principles of their disciplines to processes based on living cells or subcomponents of such cells. The problems of detailed equipment design, sensor development, control algorithms, and manufacturing strategies can utilize the principles from these disciplines. Biochemical engineering is more limited in the sense that it draws primarily from chemical engineering principles and is broader in the sense that it is not restricted to well-defined artificially constructed processes, but it can be applied to natural systems [START_REF] Dutta | Fundamentals of biochemical engineering[END_REF].

The fundamental training of biologists and engineers is distinctly different. In the development of knowledge in the life sciences, unlike chemistry and physics, mathematical theories and quantitative methods (except statistics) have played a secondary role. Most progress has been due to improvements in experimental tools. Results are qualitative and descriptive models are formulated and tested. Consequently, biologists often have incomplete backgrounds in mathematics but are very strong with respect to laboratory tools and, more importantly, with respect to the interpretation of laboratory data from complex systems. Engineers usually possess a very good background in the physical and mathematical sciences. Often a theory leads to mathematical formulations, and the validity of the theory is tested by comparing predicted responses to those in experiments. Quantitative models and approaches, even to complex systems, are strengths. Biologists are usually better at the formation of testable hypotheses, experimental design, and data interpretation for complex systems. Engineers are typically unfamiliar with the experimental techniques and strategies used by life scientists. The skills of the engineer and life scientist are complementary. To convert the promises of molecular biology into new processes, new products must require the integration of these skills. To function at this level, the engineer must have a solid understanding of biology and its experimental tools [START_REF] Shuler | Bioprocess Engineering: Basic Concepts[END_REF].

Due to the limitations of the thesis, it is not possible to cover all dimensions involved in bio-process engineering. However, we consider many topics that are used to develop the problems that are studied. In the next sections, appropriate references are suggested.

Bioreactors: Design, kinetics and operation mode

The bioreactor is a device used to carry out any kind of bioprocess; examples include fermenters or enzyme reactor [START_REF] Gold | Compendium of chemical terminology : IUPAC recommendations[END_REF]. Put simply, a bioreactor is a vessel in which a chemical process is carried out that involves organisms or biochemically active substances derived from such organisms. This process can either be aerobic or anaerobic. There are different types and sizes, depending on its use (industrial or research, for instance). Everything concern with the implementation and dynamics involved in the process is the basis for the bioreactor design.

Bioreactor design is an integral part of biotechnology, an area with rather vague and contested borders. Biotechnology is not simply the sum of microbiology, genetics, biochemistry, engineering, etc.; it is also the integration of these disciplines, and this involves quite a bit more than just simple addition. Integration and application are the keywords that can be found in most definitions of biotechnology. In particular with the design of bioreactors, integration of biological and engineering principles is essential. The bioreactor should be designed such that the specific biological and technological demands of a process are met. Naturally, the quality and price of a product are important for commercial realization. The aim of bioreactor design can thus be defined as "minimization of the costs of the pertinent product while retaining the desired quality, and this within the biological and technological constraints." This does not mean a priori that minimizing the costs of the bioreactor also means minimizing the costs of the integral process. This depends largely on the cost-determining part(s) of the process. If running the bioreactor is cost determining, then maximization of the overall volumetric productivity of the bioreactor is, in general, the rational approach. If, on the other hand, downstream processing is cost determining, then maximization of the product concentration in the bioreactor is, in general, the rational thing to do. However, here again, integration is the keyword. Bioreactor design should be an integral part of the overall process design. The words bioreactor, biocatalyst and product have been used in general terms [START_REF] Van't Riet | Basic Bioreactor Design[END_REF].

Productivity and product concentration

Overall volumetric productivity Q p [mol m -3 s -1 ] (it is also common to use a yearly basis) is the average production capacity per unit volume and time of the bioreactor. The overall volumetric productivity is confined, on the one hand, by physical constraints such as mass and heat transfer and, on the other hand, by the biocatalyst concentration x [mol m -3 ] and activity of the bio-catalyst, expressed as the substrate consumption rate -µ s [mol m -3 s -1 ]. Maximization of the overall volumetric productivity of the bioreactor in principle means minimization of the costs of investment because one can suffice with smaller equipment. It also usually means lower operating costs. In general, it also means that it is desirabled to operate the bioreactor as close as possible to the physical constraints. This physical limitation is the result of mass and heat transfer limitations, which are stoichiometrically related to product formation.

In addition to limitations by mass and heat transfer and concentration of the biocatalyst, the overall volumetric productivity of the bioreactor is determined by the overall productivity of the biocatalyst, Pr x , defined as the total moles of a product that are produced by 1 mol of biocatalyst during its operational lifetime t l [s]. Pr x is related to the specific product production rate q p [s -1 ] (moles of product produced per mol of biocatalyst per second) as

Pr x = tl 0 q p dt = - tl 0 Y ps µ s x dt. (1.1)
The definition of Y ps , the overall yield of product on substrate (total moles of product produced per total mol of substrate), leads to

Q p = 1 t l t l 0 q p xdt = 1 t l t l 0 -Y ps µ s dt. (1.
2)

The time needed to empty, clean, refill, restart, etc. the bioreactor between two operations is the so-called down-time, which is symbolized by t a [s]. In case t a is relevant, it can be introduced in the equations above by replacing 1 t l preceding the integral, by 1 t l +t a . In addition to the molar productivity used above, the mass productivity (kg product instead of mol) is also commonly used in engineering (conversion from one to the other by means of the molecular weights). It is also common practice to use hour, day or year as a unit of time. The search for and the development of a useful biocatalyst with a suitable yield, specific activity and stability is, in the beginning, the task of micro-biologists, biochemists, molecular biologists, protein engineers, etc. However, in particular with respect to stability, the process engineer also has the means available -among others immobilization-to improve the stability of biocatalysis.

The effect of the composition of the product stream leaving the bioreactor on the costs of downstream processing is substantial. Therefore, it is essential to take this into account when designing the bioreactor. In practice, this often means that the bioreactor is designed such that the concentration of the product is as high as possible. The end concentration of product C p [molm -3 ] in the bioreactor depends on µ s , Y ps , and the residence time in the bioreactor. For a batch reactor, with t f [s] as the time that the batch lasts, this leads to

C p = Y ps t f 0 -µ s dt, (1.3)
and for a continuous reactor with a liquid throughflow F l [m 3 s -1 ] and a volume V [m 3 ]

C p = -Y ps µ s V F l (1.4)
Concentration of product is a key-parameter when the downstream processing is the cost-determining part of the integral process. Product recovery is often a laborious and expensive operation, in particular when diluted aqueous solutions are involved such as are usually encountered in biotechnology. However, it has become clear that the aqueous reaction medium, which was for a long time supposed to be essential for biocatalysis, can be replaced to a large extent by a suitable organic solvent.

Volumetric productivity in a batch culture is the objective to maximize in the problem proposed and analyzed in chapter 2 of this thesis.

Enzyme kinetics

Enzymes are biological catalysts that are protein molecules in nature. They are produced by living cells (animal, plant, and microorganism) and are absolutely essential as catalysts in biochemical reactions. Almost every reaction in a cell requires the presence of a specific enzyme. A major function of enzymes in a living system is to catalyze the creation and breaking of chemical bonds. Therefore, like any other catalysts, they increase the rate of reaction without themselves undergoing permanent chemical changes.

The catalytic ability of enzymes is due to its particular protein structure. A specific chemical reaction is catalyzed at a small portion of the surface of an enzyme, which is known as the active site. Some physical and chemical interactions occur at this site to catalyze a certain chemical reaction for a certain enzyme. Enzyme reactions are different from chemical reactions, as follows [START_REF] Dutta | Fundamentals of biochemical engineering[END_REF]:

1. An enzyme catalyst is highly specific, and catalyzes only one or a small number of chemical reactions. A great variety of enzymes exist, which can catalyze a very wide range of reactions.

2. The rate of an enzyme-catalyzed reaction is usually much faster than that of the same reaction when directed by nonbiological catalysts. Only a small amount of enzyme is required to produce a desired effect.

3. The reaction conditions (temperature, pressure, pH, and so on) for the enzyme reactions are very mild.

4. Enzymes are comparatively sensitive or unstable molecules and require care in their use.

Because an enzyme is a protein whose function depends on the precise sequence of amino acids and the protein's complicated tertiary structure, large-scale chemical synthesis of enzymes is impractical if not impossible. Enzymes are usually made by microorganisms grown in a pure culture or obtained directly from plants and animals.

Enzyme kinetics adresses the rate of enzyme reaction and how it is affected by various chemical and physical conditions. Kinetic studies of enzymatic reactions provide information about the basic mechanism of the enzyme reaction and other parameters that characterize the properties of the enzyme. The rate equations developed from kinetic studies can be applied in calculating reaction time, yields, and optimum economic condition, which are important in the design of an effective bioreactor.

Assume that a substrate (S) is converted to a product (P) with the help of an enzyme (E) in a reactor as S E -→ P. If you measure the concentrations of substrate and product with respect to time, the product concentration will increase and reach a maximum value, whereas the substrate concentration will decrease. The rate of reaction µ s can be expressed in terms of either the change in substrate s or product concentrations C p as follows:

µ s = -ds dt , (1.5
)

µ s = dC p dt . (1.6)
In order to understand the effectiveness and characteristics of an enzyme reaction, it is important to know how the reaction rate is influenced by reaction conditions such as substrate, product, and enzyme concentrations. If we measure the initial reaction rate at different levels of substrate and enzyme concentrations, we obtain a series of curves and we can conclude the following:

1. The reaction rate is proportional to the substrate concentration (that is, first-order reaction) when the substrate concentration is in the low range.

2. The reaction rate does not depend on the substrate concentration when the substrate concentration is high, since the reaction rate changes gradually from first order to zero order as the substrate concentration is increased.

3. The maximum reaction rate µ max is proportional to the enzyme concentration within the range of the enzyme tested.

Considering the remarks above, the following rate equation (Henri, 1902) [START_REF] Dutta | Fundamentals of biochemical engineering[END_REF] is proposed

µ s (s) = µ max s K M +s , (1.7)
where µ max and K M are kinetic parameters that need to be experimentally determined. Eq. (1.7) expresses the three preceding observations fairly well. The rate is proportional to s (first order) for low values of s, but with higher values of s, the rate becomes constant (zero order) and equal to µ max . Since (1.7) describes the experimental results well, we need to find the kinetic mechanisms that support this equation. The reaction rate equation can be derived based on the following assumptions:

1. The total enzyme concentration stays constant during the reaction (C

E 0 = C ES + C E ).
2. The amount of an enzyme is very small compared to the amount of substrate. Therefore, the formation of the enzyme substrate complex does not significantly deplete the substrate.

3. The product concentration is so low that product inhibition may be considered negligible.

In addition to the preceding assumptions, there are three different approaches to derive the rate equation: the Michaelis-Menten approach (1913), Briggs-Haldane approach (1925) and numerical solution. We focus on the Michaelis-Menten approach, which is mostly used in this work. In the first place, we note that the Michaelis-Menten equation has the same form as the Monod equation ( 1942), but differs in that it is based on theoretical considerations, while the latter is empirical.

If the slower reaction, ES k 3

-→ P + E, determines the overall rate of reaction, the rate of product formation and substrate consumption is proportional to the concentration of the enzyme-substrate complex as

µ s = dC p dt = ds dt = k 3 C ES (1.8)
Unless otherwise specified, the concentration is expressed as molar unit, such as kmol/m 3 or mol/L. The concentration of the enzyme-substrate complex C ES , can be related to the substrate concentration s and the free-enzyme concentration C E based on the assumption that the first re-

versible reaction S + E k 1 ← → k 2
ES is in equilibrium. Then, the forward reaction is equal to

k 1 sC E = k 2 C ES .
By substituting the last equality into (1.8), the rate of reaction can be expressed as a function of s and C E , of which C E cannot be easily determined. If we assume that the total enzyme contents are conserved, the free-enzyme concentration C E can be related to the initial enzyme concentration C E 0 (see first assumption above). Thus, we now have three equations from which we can eliminate C E and C ES to express the rate expression as the function of substrate concentration and the initial enzyme concentration. By substituting the last equations for C E and rearranging for C ES , we obtain

C ES = C E 0 s k 2 k 1 +s
.

Substitution of the equation above into (1.8) results in the final rate equation

µ s = dC p dt = ds dt = k 3 C E 0 s k 2 k 1 +s = µ max s K M +s (1.9)
which is known as the Michaelis-Menten equation and is identical to the empirical expression (1.7). K M in equation (1.9) is known as the Michaelis constant.

The unit of K M is the same as s. When K M is equal to s, µ s is equal to one half of µ max according to (1.9). Therefore, the value of K M is equal to the substrate concentration when the reaction rate is half of the maximum rate µ max . K M is an important kinetic parameter because it characterizes the interaction of an enzyme with a given substrate. Another kinetic parameter in (1.9) is the maximum reaction rate µ max , which is proportional to the initial enzyme concentration. The main reason for combining two constants k 3 and C E 0 into one lumped parameter µ max is due to the difficulty of expressing the enzyme concentration in molar units. To express the enzyme concentration in molar units, we need to know the molecular weight of the enzyme and the exact amount of pure enzyme added, both of that are very difficult to determine. Since we often use enzymes which are not in pure form, the actual amount of enzyme is not known. Enzyme concentration may be expressed in mass unit instead of molar units. However, the amount of enzyme is not well quantified in mass units because the actual contents of an enzyme can differ widely depending on its purity. Therefore, it is common to express enzyme concentration as an arbitrarily defined unit based on its catalytic ability.

Michaelis-Menten kinetics (or Monod growth in the context of microalgae cultivation) is mainly used in chapters 2 and 3 of this thesis.

Batch operation mode

The simplest reactor configuration for any enzyme reaction is the batch mode. A batch enzyme reactor is normally equipped with an agitator to mix the reactant, and the pH of the reactant is maintained by employing either a buffer solution or a pH controller. An ideal batch reactor is assumed to be well mixed so that the contents are uniform in composition at all times. and

K M ln( s 0 s ) + (s 0 -s) = µ max t. (1.12)
This equation shows how s is changing with respect to time. With known values of µ max and the change in s with time in a batch reactor can be predicted from this equation.

In a batch reactor, the reactants and the catalyst are placed in the reactor, which is then closed to the transport of matter; the reaction is allowed to proceed for a given time, whereupon the mixture of unreacted material together with the products is withdrawn. A provision for mixing may be required [START_REF] Gold | Compendium of chemical terminology : IUPAC recommendations[END_REF]. A special case of batch reactor with light influence (photobioreactor) is the focus in chapter 2 of this thesis.

Continuous stirred-tank reactor

A continuous stirred-tank reactor (CSTR) is an ideal reactor that is based on the assumption that the reactor contents are well mixed. Therefore, the concentrations of the various components of the outlet stream are assumed to be the same as the concentrations of these components in the reactor. Continuous operation of the enzyme reactor can increase the productivity of the reactor significantly by eliminating the downtime. It is also easy to automate to reduce labor costs. The substrate balance of a CSTR can be set up, as follows

Fs 0 -Fs + µ s V = V ds dt , (1.13)
where F is the flow rate and V is the volume of the reactor contents. It should be noted that µ s is the rate of substrate consumption for the enzymatic reaction, while ds dt is the change in substrate concentration in the reactor. As can be seen in (1.13), µ s is equal to ds dt when F is zero, which is the case in batch operation.

For the steady-state CSTR, the substrate concentration of the reactor should be constant. Therefore, ds dt is equal to zero. If the Michaelis-Menten equation can be used for the rate of substrate consumption µ s , equation (1.13) can be rearranged as

F V = D = 1 τ = µ max s (s 0 -s)(K M +s) (1.14)
where D is known as the dilution rate, and is equal to the reciprocal of the residence time τ. It is common in biochemical engineering to use the term dilution rate rather than the term residence time, with which chemical engineers are more familiar.

Equation (1.14) can be rearranged to give the linear relationship

s = -K M + µ max sτ s 0 -s . (1.15)
Michaelis-Menten kinetic parameters can also be estimated by running a series of steady-state CSTR runs with various flow rates and plotting s versus sτ s 0 -s .

A chemostat is a CSTR for growing bacteria cells. Specifically, the chemostat is a bioreactor in which constant growth conditions for microorganisms are maintained over prolonged periods of time by supplying the reactor with a continuous input of nutrients and continuous removal of medium [START_REF] Gold | Compendium of chemical terminology : IUPAC recommendations[END_REF]. Dynamics in a chemostat are discussed in chapter 3 of this thesis.

Enzyme inhibition

Inhibitors are substances normally found in a naturally occurring reaction system, foreign-substances contaminants, or unexpected reaction retardants. As discussed previously, substrate or products may also be inhibitors. Degree of inhibition may depend on pH or the presence of other substances in the reaction mixture. The most common inhibitors that must be dealt with are reaction substrates and/or products.

Inhibition of enzymatic reactions may be reversible, irreversible, or somewhere between the two extremes. Even if inhibition is reversible, the rate of reversal may be so slow that inhibition must still be considered essentially irreversible. Either reversible or irreversible inhibition may be so slow that their effects must be considered as rate limited. Because absolutely irreversible inhibition results in permanently inactive enzymes, the process is also called deactivation. The term is not reserved for irreversible inactivation, however.

Several standard mechanisms utilized to explain the effects of reversible inhibition are based on Michaelis-Menten kinetics and can be modeled by using modifications of Equation (1.8) with linear functions of inhibitor concentration I. The modifications are as follows for the four major categories.

Competitive inhibition:

1

K app A = 1 K A 1 + 1 K IC . (1.16)
2. Uncompetitive inhibition:

1 K app 0 = 1 K 0 1 + 1 K IU .
(1.17 The overall rate equation(s) become complex when all these factors are included. Only in the simpler cases can accurate integrations be carried out without the use of numerical integration. Most initial rate equations fall into one standard category or another. The models for these categories are conveniently stated in simpler standard forms analogous to the Michaelis-Menten equation [START_REF] Mcduffie | Bioreactor Design Fundamentals[END_REF][START_REF] Shuler | Bioprocess Engineering: Basic Concepts[END_REF].

• For competitive inhibition µ s = µ max s K M 1+ I K IC +s
.

(1.18)

• For uncompetitive inhibition

µ s = µ max s K M +s 1+ I K IU . ( 1.19)
• For mixed inhibition

µ s = µ max s K M +s+ K M K IC + sI K IU . (1.20)
Where s represents the substrate concentration and I represents the inhibitor concentration. Confusing terminology has arisen for the special case of mixed inhibition, where K IC = K IU , which traditionally is called noncompetitive inhibition. Not only is there confusion between the terms noncompetitive inhibition and uncompetitive inhibition, but also among the various types of mixed inhibition, which may also have been called noncompetitive inhibition. The Nomenclature Committee of the International Union of Biochemistry recommends that this special case of mixed inhibition, where K IC = K IU , be called pure noncompetitive inhibition and that the term noncompetitive inhibition not be used at all. In line with these recommendations, the equation for pure noncompetitive inhibition is

µ s = µ max s (K M +s) 1+ I K I . (1.21)
In some cases, the maximum rate µ max may be determined by direct observation. Moreover, K M is equal to the substrate concentration giving a rate equal to half the maximum rate. Care should be taken in the use of this simplified analysis to be sure that a true maximum rate is obtained and that the reaction actually is uninhibited. Note thata the catalytic constant,K 0 , is the slope of a plot of µ max versus E 0 , for the uninhibited reaction; then

K 0 K A = K M . (1.22)
The analysis of inhibition models according to Michaelis-Menten kinetics also may be applied conveniently to cell growth inhibition because the simplified equation forms presented here are the same as those derived from the Monod equation. This result naturally follows, because the Monod equation presents a specific growth rate µ as a Michaelis-Menten function of cell substrate concentration. However, cell multiplication may require more complex models and different treatment, just as some complex enzymatic reactions do. Noncompetitive inhibition terms are used in the example studied in chapter 3 of this thesis.

Microalgae cultivation: light influence and photobioreactors

Microalgae are microscopic photosynthetic microorganisms. They are a type of eukaryotic cell, and they contain similar organelles such as chloroplasts and nuclei. Microalgae are generally more efficient than land based plants in utilizing sun light, CO2, water and other nutrients which translates into higher biomass yields and higher growth rates. They can also grow in a variety of aquatic environments and can be grown without the use of fertilizers and pesticides which results in less waste and pollution [START_REF] Richmond | Handbook of microalgal culture : biotechnology and applied phycology[END_REF].

There are many areas where microalgae are needed. The whole algal biomass can be used as a source of protein or valuable chemicals (pigments, enzymes) can be extracted. For instance, in pharmaceutical applications several studies have shown the potential of microalgae as therapeutical agents. Microalgae are already used as a source of nutritional supplements, as an additive for cosmetics, in the treatment of wastewater and as a potential source of oil for biofuels [START_REF] Bernard | Hurdles and challenges for modelling and control of microalgae for CO2 mitigation and biofuel production[END_REF].

Reduction of CO2 emissions is expedient, but even a 30% reduction (as agreed by some international legislations) is still not enough to stabilize the CO2 levels within a safe zone. This requires the development of alternative biofuels, of which biodiesel and bioethanol have the greatest market potential at the moment. Biodiesel is a popular alternative to petroleum-based diesel; it can be used in regular diesel engines. It is ecofriendly, non toxic, biodegradable and when burned, due to its low sulphur content, it produces fewer emissions than its petroleum counterpart. It can be prepared from renewable sources like edible and non edible vegetable oils, animal fats and even waste cooking oils. The most common concern with first generation biofuels is that as production capacities increase, so does their competition with agriculture for arable lands used for food production [START_REF] Castellanos | Batch and continuous studies of chlorella vulgaris in photobioreactors[END_REF].

Different approaches can be taken when looking to grow microalgae in large volume. These include outdoor ponds, with light supplied by the sun, and photo-bioreactors, which can be outdoor or indoor with light supplied by electric lights. Species control is better achieved under closed conditions that are very common in a laboratory setting. Closed systems provide a better opportunity to meet specific demands and to control and optimize cell growth parameters. In the case of photo-bioreactors, many design considerations need to be made depending on the end goal. There are basic design features that should be considered regardless of the configuration; these include materials used for the set-up, source of light, circultion of the algae through the reactor, supplying CO2, and the control of other parameters such as pH, temperature and nutrient concentration in the media [START_REF] Chisti | Biodisel from microalgae[END_REF].

To estimate growth rates, one must have a series of measurements, at different times, that will permit the calculation of the rate of change in biomass concentration. Cell number should be counted, either through a direct method, as through light microscopy with a hemocytometer, or indirectly through biomass concentration (as dry weight) or optical density as long as these measurements correlate linearly with the number of cells. Under a typical homogenous batch regime (in a closed system), microalgae (and other microorganisms) will pass through the following growth phases: The individual phases, shown in Figure 1.4, are not always clearly defined; their length or slope might change according to the culture conditions. During the lag phase (phase A) the microalgae cells are adapting to the new environment conditions, at the end of the lag phase, the cells are well adjusted and then start to multiply rapidly; this is the exponential phase (phase B. In this phase, the number of living cells, doubles regularly with time. During this period, the cell growth is described by the differential equation

dx dt = µ s x, x(0) = x 0 , (1.23)
where x is the number of cells (or cell concentration) at different times and µ s is the specific growth rate.

The deceleration growth phase follows the exponential phase. In this phase, growth decelerates due to either depletion of one or more essential nutrients or the accumulation of toxic byproducts of growth. For a typical bacterial culture, these changes occur over a very short period of time. The rapidly changing environment results in unbalanced growth. During unbalanced growth, cell composition and size will change. In the exponential phase, the cellular metabolic control system is set to achieve maximum rates of reproduction. In the deceleration phase, the stresses induced by nutrient depletion or waste accumulation cause a restructuring of the cell to increase the prospects of cellular survival in a hostile environment. These observable changes are the result of the molecular mechanisms of repression and induction. Because of the rapidity of these changes, cell physiology under conditions of nutrient limitation is more easily studied in a continuous culture.

The stationary phase (phase C) starts at the end of the deceleration phase, when the net growth rate is zero (no cell division) or when the growth rate is equal to the death rate. Although the net growth rate is zero during the stationary phase, cells are still metabolically active and produce secondary metabolites. Primary metabolites are growth-related products and secondary metabolites are nongrowth-related. In fact, the production of certain metabolites is enhanced during the stationary phase (e.g., antibiotics, some hormones) due to metabolite deregulation. The reason for termination of growth may be either exhaustion of an essential nutrient or accumulation of toxic products. If an inhibitory product is produced and accumulates in the medium, the growth rate will slow down, depending on inhibitor production, and at a certain level of inhibitor concentration, growth will stop. Ethanol production by yeast is an example of a fermentation in which the product is inhibitory to growth. Dilution of toxified medium, with the addition of an unmetabolizable chemical compound complexing with the toxin, or simultaneous removal of the toxin would alleviate the adverse effects of the toxin and yield further growth.

The death phase (phase D) follows the stationary phase. However, some cell death may start during the stationary phase, and a clear demarcation between these two phases is not always possible. Often, dead cells lyse, and intracellular nutrients released into the medium are used by the living organisms during stationary phase. At the end of the stationary phase, because of either nutrient depletion or toxic product accumulation, the death phase begins. The rate of death usually follows first-order kinetics

dx dt = -µ d x, (1.24)
where x is the concentration of cells at the end of the stationary phase and µ d is the first order death rate constant. A complete analysis of a batch growth is given in [START_REF] Shuler | Bioprocess Engineering: Basic Concepts[END_REF].

Light influence in growth and productivity

The central issue involved in mass cultivation of photoautotrophic microalgae,who are organisms that carry out photon capture to acquire energy, concerns effective use of strong light for photosynthetic productivity of cell mass and secondary metabolites. This is particularly true for mass cultivation of microalgae outdoors, in which effective use of solar energy is the foundation on which the prospects for this biotechnology rest. Light energy received by photoautotrophic microorganisms is a function of the photon flux density (PFD) reaching the culture surface. The cells absorb only a fraction of the photon flux, the actual size of which is governed by several factors, including cell density, the optical properties of the cells, length of the optical path of the reactor and rate of culture mixing. Photons that are not absorbed by the cells' photosynthetic reaction centers dissipate mostly as heat or may be reflected. As a rule, microalgal mass cultures reflect only a small or very small fraction of the photons impinging on a culture surface.

Since essentially all photons of a high flux density may be captured by high cell density cul-tures, cell density will continue to increase exponentially until all photosynthetically available photons are absorbed. Once this cell density is reached, cell mass accumulates at a constant, linear rate until light per cell or some substrate in the culture medium becomes overly low, or alternatively, some inhibitory activity or conditions arrest cell growth. In the light-limited linear growth phase, the relationship between biomass output rate and light energy absorbed by the culture (I 0 A) can be expressed as follows

I 0 A = µ s xV Y , (1.25)
where the radiation intensity at the edge of the illuminated side of the reactor is denoted by I 0 , A represent the irradiated culture area, µ s is the specific growth rate, x represent biomass density, V is the total culture volume and Y represent bioenergetic growth yield. This relationship implies that the biomass output rate in continuous cultures µ s x is determined by area-volume relationships (A/V), and in order to obtain high cell densities it is mandatory to use a reactor of high A/V ratio. Finally, if the value of Y for a particular microalga is a constant, the specific growth rate µ s can be altered by adjusting x without changing any other culture parameters.

Maximal culture productivity may be obtained only when culture nutritional requirements are satisfied and temperature is about optimal. There exists, indeed, a strong interaction between light and temperature. At higher light intensities, the photosynthetic rate increased with an increase in temperature, and at high, arround optimal temperature, the photosynthetic rate increased with increasing light intensities. Thus, an elementary aspect of the interaction of light and temperature is that the optimal temperature for photosynthesis increases with increasing light intensiy.

As a rule, growth of phototrophic mass cultures should be limited by light only. Efficient utilization of strong light by the individual cells in the culture is associated, however, with many constraints: One difficulty rests with the fact that the photosynthetic photon flux density (PPFD) required to saturate the photosynthetic units in the cell is usually 1/5 or 1/10 the PPFD impinging on the culture at midday. Even relatively short exposure of the photosynthetic unit (PSU) or photosynthetic reaction center, to a light dose much above saturation may impair the photosynthetic complex and may reduce productivity. The kinetic response of an algal cell to light intensity is shown in a generalized shape of the curve relating algal growth to the intensity of the light source (see figure 1.5), provided the light source is strictly the sole limiting factor for growth and cell development. The main features of this curve are as follows: At some very low light intensity, the resulting low growth-rate is balanced by decay and the net growth is zero (compensation point). As light becomes more intense, growth is accelerated, the initial slope of the curve representing maximal efficiency of growth in response to light. With a further increase in light intensity, the light saturation function is reached, at which point the growth rate is the maximal attainable; a further increase in light intensity above this point would not result in a further increase in the growth rate; but may become injurious, manifested initially by a decreased growth rate and culminating in photo-damage and, in extreme cases, in culture death (see Chapter 2 in [START_REF] Richmond | Handbook of microalgal culture : biotechnology and applied phycology[END_REF]). If all growth conditions are optimal and the culture is optically thin (i.e., of low cell density), the intensity of the light source is indeed the sole determinant of light availability to the cells and hence the sole factor that controls growth. Under these conditions, the effect of the intensity of the light source on algal photosynthesis and growth is faithfully portrayed by the light response curve shown in figure 1.5. It is thus essential to note at this point that much misunderstanding concerning the complex effect of light on culture growth has originated from irrelevant application of the light curve for interpretation of the growth response of mass cultures, in which cell density, as a rule, is such that mutual shading may greatly modify light availability for the individual cells. This internal shading (clearly visible in that light does not pass through the culture's optical path, being essentially fully absorbed in the outer surfaces) results in cells receiving light intermittently, a phenomenon augmented by the fact that light energy attenuates exponentially in passing through the culture column. The higher the cell density, the shorter the depth light penetrates into the culture. Two light zones are thereby established in the culture: the outer illuminated volume, in which light is sufficient to support photosynthesis (i.e., the photic zone); and the dark volume, in which net photosynthetic productivity cannot take place since light intensity is below the compensation point (figure 1.5). The higher the population density (and the longer the optical path), the more complex it becomes to address the basic requirements for efficient utilization of strong light, i.e., an even distribution of the available light to all cells in the culture, at an optimal dose per cell (to be elucidated somewhat later). Clearly then, when mutual shading prevails, cells are not exposed to continuous illumination but rather to cycles of light and darkness (L-D cycle), which may take scores of milliseconds to a few seconds to complete, depending on the optical path and the extent of turbulence in the culture. The endless combinations of light intermittency expressed in L-D cycles to which the individual cells are exposed at a given instant relate to two basic parameters: first, the ratio between the light and the dark period in the cycle and, second, the frequency of the cycle. As shall be elucidated, the higher the frequency of the L-D cycle, the more efficient strong light may be used for photosynthesis. It can be readily seen, therefore, that the effect of light on photosynthetic productivity (i.e., cell mass produced per illuminated area per time) as depicted in figure 1.5 may be misleading or altogether irrelevant; in effect, it ignores the decisive impact on productivity exerted by other factors that concern the photon flux to which the cells are exposed, i.e., cell density, the length of the optical path and the extent of culture turbulence, all of which represent major determinants affecting phototrophic productivity, not less important than the intensity of the light source. In what follows, these determinants, all having decisive effects on photosynthetic productivity of mass cultures, will be elucidated.

A complete study of light influence in microalgae growth, including growth inhibitory substances and conditions (see, for instance, the example in chapter 3), average radiation intensity (see appendix A and B in chapter 2), effective use of sunlight and high irradiance for photosynthetic productivity, photosynthetic efficiency in mass cultures and other interesting topics can be found in [84, chapter 8].

Photobioreactors

Photobioreactors (PBR) are reactors in which phototrophs (microbial, algal or plant cells) are grown or used to carry out a photobiological reaction. In a broad sense, the open shallow basins widely used for microalgae cultivation could also be viewed as photobioreactors [94].

At present, commercial production of phototrophic microbial biomass is limited to a few microalgal species that are cultivated in open ponds by means of a selective environment or a high growth rate. Most microalgae cannot be maintained long enough in outdoor open systems because of the risk of contamination by fungi, bacteria and protozoa, and competition by other microalgae that tend to dominate regardless of the original species that are used as inoculum. PBR offer a closed culture environment, which is protected from direct fallout, relatively safe from invasion by competing microorganisms, and where conditions are better controlled ensuring dominance of the desired species. Thus, PBR allow the exploitation of the potential of the more than 50 000 microalgal species known, many of which may be interesting sources of high value compounds [START_REF] Richmond | Handbook of microalgal culture : biotechnology and applied phycology[END_REF].

Photobioreactors can be defined as culture systems for phototrophs in which a great proportion of the light (> 90%) does not impinge directly on the culture surface, but has to pass through the transparent reactor's walls to reach the cultivated cells. Consequently, PBR do not allow, or strongly limit, direct exchange of gases and contaminants (dust, microorganisms, etc.) between the culture and the atmosphere.

Photobioreactors can be classified on the basis of both design and mode of operation. In design terms, the main categories of reactors are: (1) flat or tubular; (2) horizontal, inclined, vertical or spiral; and (3) manifold or serpentine. An operational classification of PBR would include (4) air or pump mixed and (5) single-phase reactors (filled with media, with gas exchange taking place in a separate gas exchanger), or two-phase reactors (in which both gas and liquid are present and continuous gas mass transfer takes place in the reactor itself). Construction materials provide additional variation and subcategories, for example (6) glass or plastic and (7) rigid or flexible PBR. Axenic PBR are reactors that are operated under sterile conditions. Although a major characteristic of PBR is their ability to limit contamination, it must be made clear that an effective barrier, and thus operation under truly sterile conditions, is not achieved except in the few special designs developed expressly for that purpose.

Design criteria for PBR should aim at achieving high volumetric productivities and high efficiency of conversion of light energy, providing, at the same time, the necessary reliability and stability to the cultivation process by means of a cost-effective culture system. An efficient PBR can not be properly designed without adequate knowledge of the physiology in mass culture of the organisms to be cultivated. Since phototrophic microorganisms are highly diverse in their morphology, nutritional and light requirements, and resistance to stresses, PBR can not be designed to adapt to all organisms and all conditions [94]. The main design criteria for PBR include: surface-to-volume ratio, orientation and inclination, mixing and degassing devices, systems for cleaning and for temperature regulation, transparency and durability of the construction material. Ease of operation and scale-up as well as low construction and operating costs also take on particular relevance in relation to commercial PBR.

There are three parameters commonly used to evaluate productivity in photobioreactors. First is volumetric productivity (VP), i.e., productivity per unit reactor volume (expressed as [g L -1 d -1 ]); the second is areal productivity (AP), i.e., productivity per unit of ground area occupied by the reactor (expressed as [g m -2 d -1 ]); the third is illuminated surface productivity (ISP), i.e., productivity per unit of reactor illuminated surface area (expressed as [g m -2 d -1 ]). It must be noted that in vertical systems the illuminated surface area comprises both the front surface receiving beam radiation and the back surface and the side walls receiving reflected and diffuse radiation. The VP is a key parameter that illustrates how efficiently the unit volume of the reactor is used. However, it should be kept in mind that VP is a function of the number of photons that enter the unit reactor volume per unit time, and as such, it is dependent on the s/v of the reactor. The higher the s/v, the higher the VP. We should also be aware of the fact that high s/v reactors may achieve high VP even if they perform poorly and that a VP of g L -1 d -1 assumes a completely different significance if obtained in a reactor of 1 cm or 5 cm light path. Care should be taken to discern between AP and ISP. In the case of ponds and horizontal or near-horizontal flat reactors, the ground surface area occupied by the system and its illuminated surface area substantially coincide, as do AP and ISP. In the case of horizontal tubular reactors, placed with tubes in contact, the illuminated surface area is 1.57 times the occupied surface area, so ISP will always be lower than AP, and both parameters can be calculated easily. In the case of horizontal tubular reactors with empty space between contiguous tubes and vertical or highly inclined systems, the situation is more complex. For example, in horizontal serpentine reactors, it is difficult to decide whether and how to compute the empty space between tubes and how to account for the fact that horizontal serpentine reactors may intercept a different proportion of the radiation impinging on the horizontal, which depends on the elevation of the sun and, hence, on the hour of the day. The performance of elevated systems may be even more difficult to evaluate, unless a fourth parameter for measuring productivity is introduced: overall areal productivity. A more complete analysis of photobioreactors is made in [START_REF] Richmond | Handbook of microalgal culture : biotechnology and applied phycology[END_REF], including scale up and comparison with open ponds.

Multiple bioreactors connected: The general gradostat

In this section, we extend the idea behind the simple chemostat to a new apparatus in order to model a property of ecological systems that is not possible to model in the simple chemostat. The idea is to capture the essentials of the new phenomenon without destroying the tractability of the chemostat either as a mathematical model or as an experimental one. Just as the chemostat is a basic model for competition in the simplest situation, the apparatus here shows promise of being a model for competition along a nutrient gradient.

The "well mixed" hypothesis for the chemostat does not allow a nutrient gradient to be generated. A basic tenet is that the nutrient concentration is the same everywhere; hence, any advantage in nutrient consumption is present everywhere. The model that incorporates a true gradient would be one involving partial differential equations; a new variable, space, must be accommodated. Systems of nonlinear partial differential equations are difficult mathematical objects to understand and analyze. Even numerical solutions pose added and significant difficulties. Moreover, even if an experimental gradient is devised, measurements that do not disturb the local environment take on new difficulties. A piece of laboratory apparatus was devised by Lovitt and Wimpenny for experiments along a nutrient gradient; this is now known as CSTR (continuous stirred-tank reactors) in series [START_REF] Mcduffie | Bioreactor Design Fundamentals[END_REF]. It is a concatenation of chemostats in which the adjacent vessels are connected in both directions. Output occurs at the first and last vessels, and those in between exchange their contents -nutrient and organisms. The flow rates in, out, and between vessels are constant and equal. The apparatus is named a gradostat. It does not occur in nature, at least in this form. Indeed, although we shall think of the apparatus as connected horizontally, the closest approximation in nature may be vertically, as in a water column. We shall see that much more imaginative connection patterns are possible. Growth along nutrient gradients does occur in abundance in nature. For example, the surface films in dental plaque represent growth along such a gradient, as does growth along the banks of a stream or along a seacoast. In a water column, sunlight replaces the nutrient as an essential source for growth.

Mathematically and experimentally there is no reason to connect the vessels linearly, to restrict the source to the left-hand vessel, or to keep the washout rates D equal so long as the volume of the fluid in each vessel is kept constant. We next describe a class of gradostat models that is sufficiently general to include all cases of biological interest and yet remain mathematically tractable. We refer to them as general gradostat.

In a loose sense the apparatus generates a "discrete" gradient; the nutrient concentrations will vary from vessel to vessel, so the "parameters" of competition change from vessel to vessel. If there is no consumption, the nutrient concentrations arrange themselves as discrete points along a linear gradient. The effect of a nutrient gradient on growth and competition can be studied with such a device. A complete stability analysis for the case of two vessels with and without competition and considering Michaelis-Menten uptake functions was explored in [START_REF] Smith | The Theory of the Chemostat[END_REF].

When the number of vessels is increased and the restriction to Michaelis-Menten uptake functions is relaxed, these computations are inconclusive. It turns out that unstable positive rest points are possible and that non-uniqueness of the coexistence rest point cannot be excluded. The main difference with the previous analysis is that we cannot exclude the possibility of more than one coexistence rest point. Let S denote the nutrient concentration and u and v the concentration of the two competitors, f u (.) and f v (.) represents growth functions of both species (Linear or Monod, for instance), y u and y v represents the yield coefficient of each species and D represents the dilution rate. Then, using the subscript i to denote concentrations of S, u, and v in vessel i = 1, ..., n, the equations take the form

             Ṡi = (S i-1 -2S i + S i+1 )D -u i y u f u (S i ) -v i y v f v (S i ) ui = (u i-1 -2u i + u i+1 )D + u i y u f u (S i ) vi = (v i-1 -2v i + v i+1 )D + v i y v f v (S i ) (1.26)
where S 0 = S(0), u 0 = v 0 = 0, S n+1 = u n+1 = v n+1 = 0, S i (0) ≥ 0, u i (0) ≥ 0, v i (0) ≥ 0, and f u , f v satisfy the following:

1. The functions f u and f v are continuously differentiable.

f u

(0) = f v (0) = 0 and f ′ u (S) > 0, f ′ v (S) > 0 for S > 0.
The Michaelis-Menten function is the prototypical example, but other functions with these properties have appeared in the literature. The principal reason for allowing quite general monotone uptake functions in this case is that we are unable to obtain sharper results under the stronger hypotheses that the uptake functions are Michaelis-Menten.

The S(0) is the input concentration of nutrient (to the leftmost vessel), and D is the dilution rate. These two parameters are under the control of the experimenter. The terms y u and y v are the yield coefficients. For convenience, one can scale substrate concentrations S i by S(0), time by 1/D (making m i , nondimensional and D = 1), and microorganism concentrations by y u S(0) and y v S(0) to obtain the less cluttered system from (1.26)

           Ṡi = S i-1 -2S i + S i+1 -u i f u (S i ) -v i f v (S i ) ui = u i-1 -2u i + u i+1 + u i f u (S i ) vi = v i-1 -2v i + v i+1 + v i f v (S i ) (1.27)
where we use the same conventions as in the unscaled equations except that S(0) = 1. Hereafter, we refer to (1.27) as the "standard" n-vessel gradostat model (see figure 1.8).

Suppose that our gradostat consists of n vessels. Let E ij be the constant (volumetric) flow rate from vessel j to i (i = j), with the convention that E ij = 0 for i = 1, ..., n. Let V i be the volume of fluid in the ith vessel, D i the flow rate from a reservoir to vessel i (D i = 0 if no such reservoir exists), S i (0) the concentration of substrate in the reservoir feeding vessel i (S i (0) = 0 if D i = 0), and C i the flow rate from vessel i to an overflow vessel (C i = 0 if no such vessel exists). The notation diag(β i ) is used to denote a diagonal matrix whose diagonal elements are given by β i ; E is the matrix of flow rates E ij .

The rate of change of the vector S(t) = (S 1 (t), ..., S n (t)) at time t in a general gradostat, in the absence of any consumers, is given by

[diag(V i )] Ṡ = ĀS + g, where Ā = E -diag[C i ] -diag( n ∑ j=1 E ji ), g = (D 1 S 1 (0), ..., D n S n (0)).
Of course, the volume V i of fluid in vessel i must be constant if this system is to describe a gradostat. This requires that

∑ j E ij + D i = ∑ i E ji + C i ,
or that the volumetric flow rates in and out of any fixed vessel be the same.

It is convenient to multiply through by [diag(V i )] -1 and obtain

Ṡ = AS + b, where A = [diag(V i )] -1 Ā and b = [diag(V i )] -1 g.
We assume that at least one vessel receives substrate, since otherwise no microorganisms can survive. Mathematically, this means that S i (0) > 0 for some i, so b = 0. From the definition of A, we have that A ii < 0 since E = 0 (excluding one trivial case of no input or output). In addition, A satisfies

A ij ≤ 0, i = j and n ∑ j=1 A ij = -V -1 i D i ≤ 0.
Our assumption that S i (0) > 0 for some i implies D i > 0, and hence strict inequality holds in the last inequality for some i.

Our principal hypothesis is that the matrix A (or, equivalently, the matrix E) will be assumed to be irreducible. Thus, the set of vessels comprising the gradostat may not be partitioned into two disjoint non-empty subsets, I and J, such that no vessel in subset J receives input from any vessel in subset I. More generally, if the matrix A does not have the property of being irreducible, one can always partition the gradostat into irreducible subsets (subgradostats) that can be studied sequentially. In this sense, there is really no loss in generality in assuming irreducibility from the start. We also mention another way to view the hypothesis of irreducibility: for any pair of distinct vessels i and j, material from vessel i can travel to vessel j, although perhaps indirectly by first passing through intermediate vessels before entering vessel j.

While we focus on irreducible gradostats, reducible gradostats may be of biological interest as well. They could be used to model a system of mountain lakes situated at different elevations, where a lake at higher elevation feeds a lake at lower elevation.

Let F u = diag[ f u (S 1 ), f u (S 2 ), ..., f u (S n ]) and let F v be defined analogously with subscript v replacing subscript u. Then, introducing consumption, the general model takes the form

   Ṡ = AS + b -F u (S)u -F v (S)v u = Au + F u (S)u v = Av + F v (S)v (1.28)
The standard model (1.26) is a special case of (1.28), where b is the vector with first component equal to 1 and all others equal to 0 and where A is the matrix with -2 in the main diagonal entries, 1 in the superdiagonal and subdiagonal entries, and 0 elsewhere. Although the standard model is of primary interest, the general gradostat described by (1.28) can be treated with the same mathematics. In figure 1.9, an irreducible gradostat configuration is described. In [START_REF] Smith | The Theory of the Chemostat[END_REF] a comprehensive analysis of the general gradostat is performed using the conservation principle and equilibrium existence and stability criteria. The gradostat is a special case of a compartmental system (see section D) and can be a motivation for the problem analyzed in chapter 4 of this thesis.

Chapter 2

Productivity optimization of microalgae cultivation in a batch photobioreactor process 2.1 Introduction

Bioreactors are laboratory and industrial devices, used for the cultivation of microorganisms. They are used in a wide variety of applications, including the production of food, beverages, pharmaceutical compounds, polymers and biofuels. In general terms, a bioreactor can be defined as a container with an inlet to introduce a cultivation media containing substances required for microorganism growth and development, and an outlet through which the produced biomass and products can be extracted. For the most part, devices are also present to enable exchange of gases such as air.

Most industrial scale bioreactors are operated under batch conditions. Thus, the bioreactor is loaded with a cultivation medium and an initial amount of biomass (inoculum), and then cells are allowed to grow for a certain period of time. Once microorganism concentration reaches a previously specified value, the reactor content is withdrawn and the unit is prepared for a new operational cycle [START_REF] Dunn | Bioreactor Modelling[END_REF]. The bioreactor characteristics depend largely on the microorganisms that need to be grown. This is clear when we conceive of photobioreactors, which are dedicated to the cultivation of microalgae [94]. A photobioreactor is a bioreactor that incorporates a source light or a disposition enabling exposure to natural solar radiation. Light is necessary to provide the input of photonic energy needed for the photosynthesis process performed by microalgae cells.

Microalgae are recognized as one of the oldest living organisms [START_REF] Falkowski | Aquatic photosynthesis[END_REF], are primitive plants (thallophytes), i.e. lacking roots, stems and leaves, have no sterile covering of cells around the reproductive cells and use chlorophyll as their primary photosynthetic pigment. While the mechanism of photosynthesis in these microorganisms is similar to that of higher plants, they are generally more efficient converters of solar energy because of their simple cellular structure. In addition, because the cells grow in aqueous suspension, they have more efficient access to water, CO2, and other substrates [START_REF] Chisti | Biodisel from microalgae[END_REF].

Microalgae can be used for the production of a wide variety of products, including high value compounds (long-chain polyunsaturated fatty acids, vitamins, and pigments) or biofuels (photobiological hydrogen, biodiesel, biomethane, bioethanol) [START_REF] Chisti | Biodisel from microalgae[END_REF]. They can also be used for environmental remediation (carbon dioxide fixation, wastewater treatment). Most of the recent interest that microalgae have received is related to their potential for biofuel production and carbon dioxide fixation based on the impact of these potential applications on the reduction of greenhouse gas emissions [START_REF] Masci | Microalgal biomass surface productivity optimization based on a photobioreactor model[END_REF]. Moreover, their high actual photosynthetic yield compared to terrestrial plants leads to large potential algal biomass production in photobioreactors of several tens of tons per hectare and per year [START_REF] Huntley | CO2 mitigation and renewable oil from photosynthetic microbes: A new appraisal[END_REF].

The photobioreactor operation is defined by the amount of light that can be provided to microalgae cells, and as biomass concentration increases, optical density of the media also increases, reducing the penetration of light in the culture an thus restricting growth. As a result, growth of microalgae decreases as time proceeds. Then, operational time in batch photobioreactors becomes a key factor determining productivity and the economic feasibility of the process. Optimization of microalgae productivity is then obviously of great interest. In recent years, many works have focused on solutions to this issue via mathematical modelling and optimal control theory [START_REF] Grognard | Optimization of a photobioreactor biomass production using natural light[END_REF][START_REF] Mairet | Modelling microalgae growth in nitrogen limited photobiorector for estimating biomass, carbohydrate and neutral lipid productivities[END_REF][START_REF] Masci | Microalgal biomass surface productivity optimization based on a photobioreactor model[END_REF]. However a more practical approach has been realized via experimental studies for specific microalgae species in different environments [START_REF] Bertucco | Continuous microalgal cultivation in a laboratory-scale photobioreactor under seasonal day-night irradiation: experiments and simulation[END_REF][START_REF] Kliphuis | Metabolic modeling of chlamydomonas reinhardtii: energy requirements for photoautotrophic growth and maintenance[END_REF][START_REF] Pruvost | Investigation of biomass and lipids production with neochloris oleoabundans in photobioreactor[END_REF].

In this work, we formulate a microalgae biomass production optimization problem in a photobioreactor in batch mode. In particular, two cases are considered separately: constant light incidence and the influence of dark/light cycles [START_REF] Janssen | Efficiency of light utilization of chlamydomonas reinhardtii under mediumduration light/dark cycles[END_REF]. To do this, we propose a simple model based on the well known Monod growth function, which is widely used in the literature [START_REF] Bernard | Hurdles and challenges for modelling and control of microalgae for CO2 mitigation and biofuel production[END_REF][START_REF] Bernard | Transient behavior of biological loop models, with application to the droop model[END_REF][START_REF] Smith | The Theory of the Chemostat[END_REF]. Mathematical results on stability and necessary optimality conditions are given for simplified nonlinear optimization problems based on the general model and its application was studied in a batch cultivation of the microalgae Chlamydomonas reinhardtii [START_REF] Harris | Chlamydomonas as a model organism[END_REF].

Optimization problem formulation 2.2.1 Model construction

In that follows, we denote by x and s the concentrations of the microorganisms and substrate respectively, Q in and Q out represent the input and output flows in the bioreactor respectively and V represent the volume of the bioreactor [START_REF] Smith | The Theory of the Chemostat[END_REF].

Assumptions 2.1

The fundamental assumptions considered in the construction of a general mathematical model of the internal dynamics in a bioreactor are:

1. The bioreactor vessel is perfectly mixed, that is, the substrate is uniformly distributed.

2. Thus, it is reasonable to assume that what is consumed c (x, s) is proportional to the amount of microorganisms, namely:

c (x, s) = µ (s) Vx, with µ (s) ≥ 0 y µ (0) = 0.
3. The growth of microorganisms is proportional to what is consumed.

4. Density of the liquid inside the reactor remain constant.

Mass balance equations for xV and sV lead us to write the general model of a bioreactor as follows:

             ds dt = Q in V (s in -s) -1 Y µ (s) x, dx dt = µ (s) x -Q in V x, dV dt = Q in -Q out , where (s, x, Q in , Q out , V) ∈ R 2 ++ × [0, Q max ] × [0, Q max ] × [0, V max ]
Q in V Input rate [h -1 ] Q out V Dilution rate [h -1 ] µ(.) Growth function [h -1 ] Y Yield constant [g cells (g substrate) -1 ]
The most common growth function is the Monod type (1942) [START_REF] Monod | La technique de la culture continue: Theorie et applications[END_REF][START_REF] Smith | The Theory of the Chemostat[END_REF] (see figure 2.1) where µ max represent the maximum specific growth rate of the microorganisms and K s represent the half saturation coefficient or Michaelis-Menten constant. For sake of simplicity, we use them in this work.

µ (s) = µ max s s + K s .

In the case where Q in = Q out = 0, reactor operated under discontinuos or batch mode. Then, the general model reduces to:

   ds dt = -1 Y µ (s) x, dx dt = µ (s) x.
(

The idea is to modify the model (2.1) from the inclusion of certain considerations about the influence of light on the dynamics, which is not considered in the system and including phenomena typical of the microalgae biochemical process.

First, we need to introduce cellular respiration. This is the sum of a series of metabolic processes that enables the microalgae to obtain the biochemical energy required for its growth. This process takes place with or without light. As a result, biomass is converted into CO2, which is represented by a term -ρx in the biomass dynamics. The parameter ρ is named respiration coefficient. Second, we consider the incident light directly influences the growth rate of microalgae biomass. This is the result of the obvious effect of light over photosynthesis, which represents the primary source of energy for microalgae. Then µ = µ(s, x, I 0 (t)) where I 0 (t) represent the time-varying incidental light (see appendix A).

From the above, the folowing model is obtained

   ds dt = -1 Y µ(s, x, I 0 (t))x, dx dt = (µ(s, x, I 0 (t)) -ρ) x.
(2.2)

Throughout this paper we consider the following assumptions (see Appendix B and C for details): Assumption 2. 2 The influence of incidental light I 0 (t) is given by

I 0 (t) =    I 0 if t 2k ≤ t < t 2k+1 ,
(light phase),

0 if t 2k+1 ≤ t < t 2k+2 , (dark phase),
with k = 0, 1, .... This mathematical simplification is enough to analize the effect of the light in the biomass production.

Assumption 2.3

The substrate input concentration should always be kept very large so as to always keep the substrate in the region where s/(s

+ K s ) ≈ 1.
That is, we do not consider substrate limitation during the studied process.

Consequently, growth function µ can be simplified by another similar increasing bounded function µ(x, t) = µx(t) c+x(t) (considering the obscuration effect over the specific growth rate) and the system (2.2) is reduced to

dx dt = µ(t)x c+x -ρx, ( 2.3)
defined by a first order nonlinear differential equation with discontinuous righthand side [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides: Control Systems[END_REF],

where the function µ (t) can be written as

µ (t) =    µ if t 2k ≤ t < t 2k+1 , (light phase), 0 if t 2k+1 ≤ t < t 2k+2 , (dark phase), (2.4)
with k = 0, 1, .... In terms of the problem studied, µ (t) = µ during the light phase (day) and µ (t) = 0 at the dark phase (night). Here t 0 > 0 denote the initial time of the process (See apendix B for details).

Problem statement

The net rate of production determines how much product (volumetric biomass in our case) one can obtained per unit time. So, in the optimization problem we can use the net rate of production as objective function. The net rate of production for a batch reactor is the quantity of product generated per batch divided by the sum of the final batch process time T and the turnaround time t a , in case of batch systems, turnaround time will include time taken in forming batches, batch execution, cleaning and forming a new batch.

The goal then would be to maximize the net rate of production in the photobioreactor in the final batch process time subject to (2.3) considering two different light environments: constant light and dark/light cycles. We can define the objective function to maximize

J (x 0 , T) = x(T)-x 0 T+t a . (2.5)
In the equation (2.5), J(x 0 , T) is the net rate of biomass production, x(T) = x(x 0 , T) represents the biomass concentration at the end of the batch, x 0 ≥ 0 represents the biomass concentration at the start of the batch process, T ≥ 0 is the final batch process time and t a > 0 is the turnaround time. The function (2.5) is called mean biomass volumetric productivity [START_REF] Pruvost | Investigation of biomass and lipids production with neochloris oleoabundans in photobioreactor[END_REF].

The reactor operator must specify the operating schedule for the operation of the reactor, and this sets the batch processing time. Doing so involves a trade-off that must be considered during the design of a batch reactor, and particularly during the specification of the operating schedule. Longer operation time implies reduction of productivity due to the optical density of the culture, so biomass reduces its growth when time passes by light obscuration effect. For this reason an optimization approach is necessary.

The dynamics described by the first order differential equation in (2.3) with the objective function defined in (2.5) are used to define the following optimization problem max

(x 0 ,T) ∈ R 2 + J (x 0 , T) = x(x 0 ,T)-x 0 T+t a , ( 2.6)
where x(x 0 , t) represents the solution of the discontinuous differential equation (2.3) with initial condition x(t 0 ) = x 0 .

Main results

In this section we introduce the first results about the problems defined in the previous part, in order to make a comparison between the optimal values of each case.

Constant light

In first place, some results will be shown for the ordinary differential equation

dx dt = µx c+x -ρx, (2.7) with initial condition x(t 0 ) = x 0 ∈ R + .
Remark Since f (x, t) = µx c+xρx is continuous as a function respect to x ∈ R, for each (x 0 , t 0 ) ∈ R 2 let x(x 0 , t) a solution of the differential equation (2.7) with initial condition x(t 0 ) = x 0 , we have that x(x 0 , t) has continuous partial derivatives with respect to x 0 [START_REF] Perko | Differential Equations and Dynamical Systems[END_REF] and ∂x ∂x 0 (x 0 , t) is the solution of the problem (variational equation) [START_REF] Hirsch | Differential Equations, Dynamical Systems, and an Introduction to Chaos[END_REF]

ẏ = µc (c+x) 2 -ρ y, y(t 0 ) = y 0 .
For simplicity, onwards we denote x(t) a solution of equation (2.7) with x(t 0 ) = x 0 , however we will return to the initial condition dependence notation if needed.

Lemma 2.4

In the nonlinear differential equation (2.7) the following holds a) The differential equation has a unique solution in the interval [0, T] with T ∈ R + 0 and this solution is bounded. b) If μ ≤ ρc then the origin is the only equilibrium point and is asymptotically stable for all x 0 > 0.

Otherwise, if μ > ρc then the origin is unstable and there is a positive and stable equilibrium point given by x e = μ ρ -c. c) A trajectory x(t) at the time t ∈ [0, T] is obtained from the implicit equation c ln

x(t) x 0 -μ ρ ln μ-ρ(c+x(t)) μ-ρ(c+x 0 ) μ-cρ = t. (2.8)
PROOF. The proof of a) follows directly from Picard-Lindelöf theorem, where the function f (x) = µx c+xρx is Lipschitz continuous with respect the variable x. Its clear also that M = max{x 0 , μ ρ -c} is an upper bound for this solution.

For prove b) in first place we note that the equilibrium points can be obtained by solving the equation µx c+xρx = 0, whose solutions are x 0 = 0 and x e = μ ρ -c. Now, to analyze its stability, we will consider the behavior of solutions for different parameter values. Suppose μ ≤ ρc, then it follows that at all x ≥ 0 µx c+xρx ≤ 0, with this dx dt ≤ 0, in this case the solution of the differential equation is decreasing in time. Also, the origin is a equilibrium solution and the uniqueness result above implies lim t→∞ x(t) = 0. Since we consider only solutions with initial conditions x (t 0 ) = x 0 > 0, it is clear from the above that if μ ≤ ρc the origin is a stable equilibrium point for all initial condition. Suppose now μ > ρc, we have that if

x 0 > µ ρ -c then µx 0 c+x 0 -ρx 0 < 0,
and

dx(t 0 ) dt < 0.
On the other hand, if x 0 < µ ρ -c using a similar analysis, we have dx(t 0 ) dt > 0. Then, if µ > ρc the origin is unstable and the equilibrium point x e is stable for any positive initial condition. This result is shown graphically in the phase plane of Figure 2.2. Finally, c) follows directly from integration, where (2.7) is separable.

In this case J(x 0 , t) is continuous and differentiable. Suppose that μ ≤ ρc, then x(x 0 , t) -x 0 ≤ 0, ∀t ≥ 0, i.e., J(x 0 , t) ≤ 0, ∀t ≥ 0 and the only possible solution of the optimization problem is the trivial (x 0 , T) = (0, 0). So, a necessary condition for the existence of a non-trivial solution of the optimization problem (2.6) is μ > ρc and this implies lim t→∞ x(t) = x e , i.e., must exist a positive and stable equilibrium point (which is ensured by Lemma 1).

On the other hand, we can note that lim t→∞ J(x 0 , t) = 0, i.e., when there are very large time values, the mean volumetric productivity is dissipated, by foregoing high mean volumetric productivity can not be achieved at very long periods of time. This implies that the optimal final time of the batch process can not be very prolonged. Therefore, defining the set

Ω(τ) = {(x 0 , t) ∈ R + 0 x 0 ≤ x e , t ≤ τ} ,
for each τ > 0 this set is compact and the Weierstrass extreme value theorem assure the existence of a solution of the optimization problem on this set. From the above, we can find the solution of the optimization problem through the deductive method [START_REF] Clarke | Functional Analysis, Calculus of Variations and Optimal Control[END_REF].

Proposition 2.5 Necessary and sufficient conditions for local optimality a) (First order) For the optimization problem defined in (2.6) considering constant light, if the pair

(x * 0 , T * ) is an optimal solution then ∂x ∂x 0 (x * 0 , T * 0 ) = 1 and J(x * 0 , T * ) = f (x(x * 0 , T * ), T * ), where f (x, t) = µx c+x -ρx. b) (Second order) Assuming a), the conditions ∂ 2 x ∂x 2 0 (x * 0 , T * ) < 0 and ∂ 2 x ∂T 2 (x * 0 , T * ) < 0,
are sufficient to ensure that a maximum is reached at that point.

PROOF. a) The result follows directly from Fermat rule [START_REF] Clarke | Functional Analysis, Calculus of Variations and Optimal Control[END_REF], where (x * 0 , T * ) must be a critical point of the objective function and therefore must fulfill

∇ x 0 ,T J(x * 0 , T * ) = 0. In particular ∇ x 0 ,T J(x 0 , T) =    ∂J ∂x 0 ∂J ∂T    .
Now, from the following equalities (stationary conditions)

∂J ∂x 0 = 1 T+t a ∂x(x 0 ,T) ∂x 0 -1 = 0, ∂J ∂T = 1 T+t a ( f (x(x 0 , T), T) -J(x 0 , T)) = 0, statement values are obtained.
b) The idea is to use the Hessian condition H(x 0 , T) for a maximization problem, which requires that the matrix H(x * 0 , T * ) be negative semidefinite in the critical point (x * 0 , T * ). Evaluating the Hessian in the critical point and assuming that the conditions in part a) are fulfilled

H(x * 0 , T * ) =    1 T * +t a ∂ 2 x ∂x 2 0 (x * 0 , T * ) 0 0 1 T * +t a ∂ 2 x ∂T 2 (x * 0 , T * )    .
If the conditions

∂ 2 x ∂x 2 0 (x * 0 , T * ) < 0, ∂ 2 x
∂T 2 (x * 0 , T * ) < 0, are satisfied, the Hessian matrix H(x * 0 , T * ) is negative semidefinite and the critical point (x * 0 , T * ) satisfies the second order necessary condition for a local maximum of the optimization problem.

As shown in Proposition above, to find it optimal candidates is necessary to solve nonlinear inequalities. However, it is not possible to do explicitly, so it is necessary to implement a numerical algorithm that allows finding an approximation to the problem solution.

Dark/light cycles

We now consider some theoretical aspects about the solutions of a nonlinear differential equation with discontinuous righthand side [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides: Control Systems[END_REF]

dx dt = µ(t)x c+x -ρx, ( 2.9)
where

µ (t) = µ if t 2k ≤ t < t 2k+1 , (light phase) 0 if t 2k+1 ≤ t < t 2k+2 , (dark phase) (2.10)
for k = 0, 1, ..., with initial condition x(0) = x 0 , which is shown in (2.3) (see appendix B).

Lemma 2.6 Carathéodory solutions

a) The differential equation with discontinuous righthand side (2.9) has unique solution in the Caratheodory sense on the interval [0, T] with T ∈ R + 0 .

b) The solutions of (2.9) in the interval [0, τ] with τ > 0, have the form

x(t) = x(t) if t 2k ≤ t < t 2k+1 , x(t 2k+1 )e -ρ(t-t 2k+1 ) if t 2k+1 ≤ t ≤ t 2k+2 ,
where

0 = t 0 < t 1 < ... < t 2k < t 2k+1 < t 2k+2 < ... < τ and x(t) is implicitly determined in t ∈ [t 2k , t 2k+1 ] (k = 0, 1, ...) by c ln x(t) x(t 2k ) -μ ρ ln μ-ρ(c+ x(t)) μ-ρ(c+x(t 2k )) μ-cρ = t -t 2k .
PROOF. a) Suppose there exists b > 0 such that the function

f (t, x) = µ(t)x c+x -ρx is defined on R = {(t, x)/t ∈ [0, a], | x -x 0 |≤ b}. This is possible, where x = x(t) is bounded on the interval [0, T].
We will prove that the function f = f (t, x) satisfies the three Carathéodory's conditions. First, f must be continuous with respect to x for almost everywhere (onwards we will denote a. e.) t ∈ R. Let t = t * fixed, the function f (x) = µ(t * ) x c+xρx is continuous throughout its domain of definition.

Secondly, the function f is measurable with respect to t for each x ∈ R. If t 2k ≤ t < t 2k+1 then f (t) = μ x c+xρx which is continuous and if t 2k+1 ≤ t < t 2k+2 then f (t) = -ρx -Dx which is also continuous, i.e. there be only discontinuities in t 2k+i with k = 0, 1, ... and i = 1, 2, 3, 4. f is therefore continuous a.e. in t and [0, T] is a Borel set, then f is Lebesgue measurable respect to t.

Now, it must be | f (t, x) | = | µ(t) x c+x -ρx | = µ(t) c+x -ρ x ≤ | µ(t) c -ρ || x | = | µ(t) c -ρ || x -x 0 + x 0 | ≤ | µ(t) c -ρ | (b + x 0 )
then by defining m(t) =|

µ(t) c -ρ | (b + x 0 ) the function m(t) is summable and | f (t, x) |≤ m(t)
, therefore f satisfies the Carathéodory's conditions in R which ensures the existence of the solution.

To proof the uniqueness it is sufficient to prove the generalized Lipschitz condition. Let (t, x), (t,

y) ∈ R | f (t, x) -f (t, y) | = cµ(t) (c+x)(c+y) -ρ (x -y) ≤ µ(t) c -ρ (x -y) ≤ µ(t) c -ρ |x -y| namely l(t) = µ(t)
cρ , then l is summable and this ensures that the solution of the Carathéodory equation (2.9) exists and is unique.

b) The first part of the proof follows directly from Lemma 2.4 part c), considering the interval [t 2k , t 2k+1 ]. The continuity of the solutions at the points t 2k+1 (k = 0, 1, ...) is ensured by Carathéodory conditions proved in Lemma 2.6. In particular, any solution is absolutely continuous and satisfies the differential equation except for a measure-zero set.

Now, if t ∈ [t 2k+1 , t 2k+2 ] the equation (2.9) is reduced to dx dt = -ρx, x(t 2k+1 ) = x(t 2k+1 ),
which can be solved by separation of variables.

Remark

1. When x(t) is a continuous piecewise-defined function, for k ∈ N fixed, is easy to show the existence of lateral derivatives

D -x(t 2k+1 ) = lim t→t 2k+1 - x(t)-x(t 2k+1 ) t-t 2k+1 = μx(t 2k+1 )-ρx(t 2k+1 )(c+x(t 2k+1 )) c+x(t 2k+1 ) , D + x(t 2k+1 ) = lim t→t 2k+1 + x(t)-x(t 2k+1 ) t-t 2k+1 = -ρx(t 2k+1 ),
analogously, for t 2k we obtain a similar result.

2. We use the definition of stability given by [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides: Control Systems[END_REF]. A solution x = φ(t) of a differential equation ẋ = f (t, x), is called stable if for each ǫ > 0 there exists δ > 0 which possesses the following property. For each x0 such that | x0φ(t 0 )| < δ, each solution x(t) with the initial data x(t 0 ) = x0 for t 0 ≤ t < ∞ exists and satisfies the inequality

| x(t) -φ(t)| < ǫ (t 0 ≤ t < ∞). 3. The point x = p is called stationary if it is a trajectory, that is, if x(t) = p is a solution of the differential equation ẋ = f (t, x)
. The term "singular point" is not used here since beside stationary points we also consider some other singular points, for instance, branching and joining points of trajectories [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides: Control Systems[END_REF]. We consider the set M = [U, L] and dist(y, M) = in f {d(y, m)/m ∈ M} the distance from y to the set M, then dist( x(t 2k+1 ), M) → 0, dist( x(t 2k+2 ), M) → 0 when k → ∞ and this implies that ∀ǫ > 0, ∃ δ > 0 such that for t 0 ≤ t < ∞ each solution x(t) with initial condition x(t 0 ) = x 0 from the δ-neighbourhood of the set M exists and satisfies the inequality dist(x(t), M) < ǫ and the set M is stable. The oscillatory behaviour follows directly from the definition of the function µ(t). b) We will use a stability result from direct Lyapunov method for discontinuous differential equations [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides: Control Systems[END_REF]. Let us consider the function

Lemma
v(t, x) = ln(1 + x),
then, in order to check the fulfillment of the stability theorem conditions, it enough guarantee that

dv dt ≡ v t + ∇v • f ≤ 0, |x| < ε 0 , (ε 0 > 0)
only in the domains of continuity of the function f (t, x).

Clearly v(t, x) ∈ C 1 , v(t, 0) = 0 and by the definition of the function, v(t, x) > 0, ∀x = 0. Now

v t + ∇v • f = ∂v(t,x) ∂x • dx dt = 1 1+x • dx dt = x 1+x µ(t) c+x -ρ < x 1+x μ c -ρ .
From the hypothesis, we have that

∂v(t,x) ∂t < x 1+x μ c -ρ < 0, ∀x = 0,
thus the function v(t, x) satisfies the condition and the origin is an stable stationary point. The Clarke generalized Jacobian [START_REF] Clarke | Functional Analysis, Calculus of Variations and Optimal Control[END_REF] (gradient in our case) of a function g is defined as the convex hull of the B-subdifferential ∂ B g, i.e., the outer limits of ∇g(x, y) when (x i , y i ) → (x, y) (The meaning is the following: consider any sequence (x i , y i ) converging to (x, y) while avoiding points at which g is not differentiable, and such that the sequence ∇g converges; then the convex hull of all such limit points is ∂g). It is possible to characterize the Clarke generalized gradient in terms of the generalized directional derivatives J 0 ((x 0 , T); (v 1 , v 2 )), in this case

∂J(x 0 , T) = (x 1 , x 2 ) ∈ R 2 : J 0 ((x 0 , T); (v 1 , v 2 )) ≥ (x 1 , x 2 ), (v 1 , v 2 ) , ∀ (v 1 , v 2 ) ∈ R 2 .
Proposition 2.8 Properties of the objective function The function J(x 0 , t) = x(t)-x 0 t+t a where x(t) is the solution of (2.9) with x(0) = x 0 is continuous, locally Lipschitz and has generalized directional derivatives in their whole domain.

PROOF. The continuity is direct from the definition and the continuity of the solution x(t). Let us now see that it is locally Lipschitz. In the first place, x(t) is locally Lipschitz with respect to t (where f (t, x) defined in the proof of lemma 2.6 is bounded). Now consider (x 0 , t), ( x0 , t0 ) ∈ R 2 , then

|J(x 0 , t) -J( x0 , t0 )| = | x(t)-x 0 t+t a -x( t)-x 0 t+t a | = |(x(t)-x 0 )( t+t a )-(x( t)-x 0 )(t+t a )| |t+t a || t+t a | = |x(t) t-x( t)t+ x 0 t-x 0 t+(x(t)-x( t))t a +( x 0 -x 0 )t a | |t+t a || t+t a | = |(x(t)-x( t)) t+( t-t)x( t)+( x 0 -x 0 )t+(t-t)x 0 +(x(t)-x( t))t a +( x 0 -x 0 )t a | |t+t a || t+t a | ≤ 1 |t+t a || t+t a | ((| t| + |x( t)| + |x 0 | + |t a |)|t -t| + (|t| + |t a |)|x 0 -x0 |) ≤ K x 0 ,t, x0 , t (x 0 , t) -( x0 , t)
where

K x 0 ,t, x0 , t = max{| t|+|x( t)|+|x 0 |+|t a |,|t|+|t a |} |t+t a || t+t a | .
Finally, to prove that J has directional derivatives, let (x 0 , t) ∈ R 2 , then

J 0 ((x 0 , t); (v 1 , v 2 )) = lim ( x0 , t)→(x 0 ,t) s→0 J(( x0 , t)+s(v 1 ,v 2 ))-J( x0 , t) s = lim ( x0 , t)→(x 0 ,t) s→0
1 s x( x0 +sv 1 , t+sv 2 )-( x0 +sv 1 ) t+sv 2 +t a -x( x0 , t)-x0 t+t a = lim ( x0 , t)→(x 0 ,t) s→0
1 s (x( x0 +sv 1 , t+sv 2 )-x( x0 , t))( t+t a )-s(( t+t a )v 1 +(x( t)-x0 )v 2 )) ( t+sv 2 +t a )( t+t a ) = lim ( x0 , t)→(x 0 ,t) ∂x ∂t ( t+t a )v 2 + ∂x ∂x 0 ( t+t a )v 1 -( t+t a )v 1 -(x( t)-x0 )v 2 ( t+t a ) 2 = 1 t+t a ∂x ∂x 0 -1, ∂x ∂t -J(x 0 , t) , (v 1 , v 2 )
Now, where

x(t) is differentiable ∀t ∈]t 2k , t 2k+1 [∪]t 2k+1 , t 2k [ with k = 0, 1, 2, ..., then J is differen- tiable ∀(x 0 , t) ∈ R×]t 2k , t 2k+1 [∪]t 2k+1 , t 2k [ and J 0 ((x 0 , t); (v 1 , v 2 )) = ∇J in this sets. Suppose that (x 0 , t 2k+1 ) ∈ R × {t 2k+1 } for k ∈ N fixed.
In these cases we have

∂J(x 0 , t 2k+1 ) = co          1 t 2k+1 +t a      dx(x 0 ,t 2k+1 ) dx 0 -1 µx(t 2k+1 )-ρ(c+x(t 2k+1 )) c+x(t 2k+1 ) -J(x 0 ,t 2k+1 )      , 1 t 2k+1 +t a      dx(x 0 ,t 2k+1 ) dx 0 -1 -ρx(t 2k+1 )-J(x 0 ,t 2k+1 )               = 1 t 2k+1 +t a dx(x 0 ,t 2k+1 ) dx 0 -1 y : y ∈ [a, b] ⊂ R ,
(co represent the convex hull) where

a = -ρx(t 2k+1 ) -J(x 0 , t 2k+1 ), b = μx(t 2k+1 )-ρ(c+x(t 2k+1 )) c+x(t 2k+1 ) -J(x 0 , t 2k+1 ).
A similar result is obtained for (x 0 , t 2k ) ∈ R × {t 2k } for k ∈ N fixed.

Remark From the result above, Rademacher's theorem [START_REF] Clarke | Nonsmooth Analysis and Control Theory[END_REF] ensures that J(x 0 , t) is differentiable almost everywhere in R 2 (the set where the function is not differentiable form a set of Lebesgue measure zero). It is ∂J = {∇J} almost everywhere in R 2 , where ∂J represent the Clarke generalized gradient of J [START_REF] Clarke | Functional Analysis, Calculus of Variations and Optimal Control[END_REF].

The problem (2.6) together with the solution of (2.9)-(2.10) is a nonsmooth optimization problem. Recall that a critical point (x * 0 , T * ) of the problem (2.6) must satisfy

∂J(x * 0 , T * ) ∋ 0, (2.11)
where ∂J(x 0 , T) represent the Clarke generalized gradient [START_REF] Clarke | Functional Analysis, Calculus of Variations and Optimal Control[END_REF][START_REF] Mordukhovich | Variational Analysis and Generalized Differentiation I: Basic Theory[END_REF].

Proposition 2.9 First-order necessary conditions An optimal solution (x * 0 , T * ) of (2.6) satisfies T * = t * 2k+1 for some k ∈ N and following conditions i) ∂x ∂x 0 (x * 0 , t * 2k+1 ) = 1, and ii) 0 ∈ 1

t * 2k+1 +t a [a * , b * ] where a * = -ρx(x * 0 , t * 2k+1 ) -J(x * 0 , t * 2k+1 ) and b * = μx(x * 0 ,t * 2k+1 )-ρ(c+x(x * 0 ,t * 2k+1 )) c+x(x * 0 ,t * 2k+1 ) -J(x * 0 , t * 2k+1 ).
It is noted that the batch process should always end at the end of a period of light.

PROOF. If the optimal solution (x * 0 , T * ) is attained in the points where the function is differentiable in the classical sense, when we evaluate the gradient in this point it should be equal to zero (Euler's necessary condition), which does not occur, i.e., the candidates to optimal solutions are the critical points (x * 0 , t * 2k+1 ) ∈ R × {t 2k+1 } for k = 0, 1, ..., such that (0, 0) t ∈ ∂J(x * 0 , t * 2k+1 ) and from the proposition above

∂J(x * 0 , t * 2k+1 ) = 1 t * 2k+1 +t a dx(x * 0 ,t * 2k+1 ) dx 0 -1 y : y ∈ [a, b] ⊂ R where a * = -ρx(x * 0 , t * 2k+1 ) -J(x * 0 , t * 2k+1 ) and b * = μx(x * 0 ,t * 2k+1 )-ρ(c+x(x * 0 ,t * 2k+1 )) c+x(x * 0 ,t * 2k+1 ) -J(x * 0 , t * 2k+1 ), then, if a critical point (x * 0 , t * 2k+1 )
is optimum, then it must fullfil the conditions in the statement.

Where the objective is a non-differentiable function, it is not possible to apply the classical direct gradient-based methods for resolution of optimization problems (which are quite efficient in general). For this reason, in the next section we will deal with an approach based on domain discretization.

Numerical approach

Throughout this section we fix t 0 = 0. In order to solve numerically problem (2.6) we consider first a tournaround time t a = 1 (in days) and define a grid of values x 0 ∈ [0, x0 ] , T ∈ [0, T], with x0 = T = 20, and step-size h = 0.1. For each of these pairs, we obtain x(T) by solving the implicit equation (2.8), and then we compute J(x 0 , T). This leads to the construction of a surface formed by the triplets (x 0 , T, J (x 0 , T)) for which the maximum of values J (x 0 , T) is found by inspection.

Chlamydomonas reinhardtii study case

Chlamydomonas reinhardtii is a single cell green alga, belonging to the chlorophytes, a group of highly adaptable species that lives in many different environments throughout the world. C. reinhardtii usually derives energy from photosynthesis, but thrives in total darkness when provided with an alternative carbon source. C. reinhardtii has been studied extensively in the past decades. It is regarded as a model organism for green microalgae because of its diverse metabolism and its ability to grow photoautotrophically as well as heterotrophically on acetate. In addition, C. reinhardtii is able to accumulate starch and produce hydrogen when grown anaerobically [START_REF] Kliphuis | Metabolic modeling of chlamydomonas reinhardtii: energy requirements for photoautotrophic growth and maintenance[END_REF]. a) (Constant light) For the algae C. reinhardtii we have the following parameter values (see ap-

pendix C) µ = 2.34 [d -1 ], ρ = 0.34 [d -1 ], c = 0.253 [g.L -1 ]
. For these values, the differential equation (2.7) has a stable positive equilibrium point in 6.63 [g.L -1 ] approx. (Lemma 1). Solving (using ❢♠✐♥❝♦♥ in ♠❛t•❛❜ for example or an interior point algorithm), the objective function reaches its maximum value 0.9389 [g.L -1 .d -1 ] and is reached in the initial condition x * 0 = 0.2 [g.L -1 ], the terminal time is T * = 2.7 [d], i.e., in day 3 of the process, the net rate of production is maximized in the case of light constant, starting with a relatively low amount of biomass and the final concentration is x(T * ) = 3.674 [g.L -1 ] approx. The pair (x * 0 , T * ) = (0.2, 2.7) satisfies the necessary and sufficient optimality condition (proposition 2.5). b) (Dark/light cycles, summer period) For the same parameter values, but considering the equation (2.9) in intervals of 14 hours of light (0.6 [d] approx.) and 10 hours of dark (0.4 [d] approx.), the objective function reaches its maximum value 0.6188 [g.L -1 .d -1 ] and is reached in the initial condition x 0 = 0.3 [g.L -1 ] and in the terminal time T = 1.6 [d], i.e., between the first and the second day of the process, the net rate of production is maximized in the case of dark/light cycles, starting with a relatively low amount of biomass, but more than in the previous case and the final concentration is x(T * ) = 1.91 [g.L -1 ] approx. Since 0 ∈ [-1.115, 1.164] then the optimal values satisfies the first order optimality condition (proposition 2.9). c) (Dark/light cycles, regular interval) For the same parameter values, but considering the model (2.9) in regular intervals of 12 hours (0.5 [d]), the objective function reaches its maximum value 0.5240 [g.L -1 .d -1 ] and is reached in the initial condition x 0 = 0.4 [g.L -1 ] and in the terminal time T = 1.5 [d], i.e., between the first and the second day of the process, the net rate of production is maximized in the case of dark/light cycles, starting with exactly the same amount of biomass that in the previous case and the final concentration is x(T * ) = 1.7099 [g.L -1 ] approx. Since 0 ∈ [-0.506, 0.441] then the optimal values satisfies the first order optimality condition (proposition 2.9). ] and is reached in the initial condition x 0 = 0.8 [g.L -1 ] and in the terminal time T = 0.4 [d], i.e., the first day of the process (during the day, in presence of light), the net rate of production is maximized in the case of dark/light cycles, starting with a relatively high amount of biomass compared to the previous example and the final concentration is x(T * ) = 1.4083 [g.L -1 ] approx. Since 0 ∈ [-0.913, 1.21] then the optimal values satisfies the first order optimality condition (proposition 2.9). turnaround time Now we will consider optimal conditions through varying incidental light, the model parameters are different respect the values used in the above simulations (but considering the same measurement units).

t a = 1 Light environment x * 0 [gL -1 ] T * [d] x(T * ) [gL -1 ] J * [gL -1 d -1 ]
For the parameters values in table 2.3, we obtain the optimal values in constant light condition using ❢♠✐♥❝♦♥ of ▼❛t•❛❜ in the optimization problem (2.6).

Discussion

In this chapter was formulated a system based on a batch photobioreactor model. From this model a nonlinear optimization problem was defined and necessary conditions for maximizing the biomass surface productivity was obtained, subject to the choice of initial biomass microalgae concentration and final operation time. In order to get analytical results we try to keep a balance between model simplicity (to handle mathematical analysis) and model complexity to capture the main phenomena that influence a photobioreactor's productivity.

We have shown numerically that, because of the day-night behaviour, the productivity rate cannot be as high as it could have been without it (constant light case). However, when the maximal growth rate is sufficiently larger than the respiration rate (for example, in summer period during light phase), we can obtain conditions for which the productivity rate is relatively close to Three practical points would be nice to highlight, which have been determined from the mathematical analysis in section 2.3:

1. There is a condition that maximizes productivity, mainly due to the obscuration effect that makes the growth rate decreases with increasing biomass concentration.

2. It is mentioned here that more light (higher intensity or longer period lighting) improve productivity, which is obvious in practice.

3. When operating in natural photoperiods, there is a moment in the day that should harvest the product of the reactor. This is the end of the day, because during the dark period, respiration causes biomass concentration decrease. The result of the optimization problem coincides with this fact.

In the particular case of microalgae Chlamydomonas reinhardtii estimates of growth parameters and respiration were obtained, the model was applied and the results were consistent with the previous analysis (section 2.3). It is interesting that during dark/light cycles in regular intervals the mean volumetric productivity decreases, being slightly more than half over a 12-hour light photoperiod compared to the case of constant light and in the same way, the terminal batch time also reduces almost by half in 12-hour photoperiod. Then (in theory) in the same time period in which mean volumetric productivity is maximized in constant light environment, two batch processes can be performed at regular intervals (12-hour photoperiod) where, under optimal conditions, higher mean volumetric productivity would be obtained in the addition of both batch process, however, this policy could be not profitable in practice. By varying the incidental light (and thus the model parameters) small variations are observed in optimal decision variables (x * 0 , T * ) in constant light case, but there exists a difference in the final The figure shows the results of the table 2.4. It may be noted that the higher the incidencial light, higher the final biomass concentration, while the terminal batch time decreases and the mean volumetric productivity increases, the initial concentration needed to achieve optimal productivity decreases slightly. concentrations and optimal mean volumetric productivity. We can see that the mean volumetric productivity in summer period (14 light hours) and at regular intervals (12 light hours) for incidental light 1000 [µmol. m -2 .s -1 ] is similar to that obtained in a batch process with constant light with incidental light 400 -500 [µmol. m -2 .s -1 ]. Are also similar the final concentrations in these cases, however, the time required to reach this level of productivity would be lower in photoperiods.

Appendices A Considering light incidence

We represent light attenuation following an exponential Beer-Lambert law

I (xz) = I 0 e -axz , (2.12)
where the attenuation at some depth z comes from the total biomass xz per surface unit contained in the layer of depth [0, z], I 0 represent the incident light and a is a light attenuation coefficient.

In microalgae, chlorophyll is mostly the cause of this shadow effect and, in model (2.1), it is best represented by a fixed portion of the biomass [START_REF] Bernard | Transient behavior of biological loop models, with application to the droop model[END_REF].

Finally, about the light source variation, will be introduced the incident light I 0 = I 0 (t) like time-varying. With such an hypothesis on the light intensity that reaches depth z, growth rates vary with depth: in the upper part of the reactor, higher light causes higher growth than in the bottom part.

Following the methodology used in [START_REF] Masci | Microalgal biomass surface productivity optimization based on a photobioreactor model[END_REF] for the simplification of the influence of light on the dynamics, supposing that light attenuation directly affects the maximum growth rate [START_REF] Huismann | Principles of the light-limited chemostat: theory and ecological applications[END_REF], the growth rate for a given depth z can then be written as

µ z (s, I (xz, t)) = µI(xz,t) I(xz,t)+K I s s+K s , with I (xz, t) = I 0 (t) e -axz .
Based on the above, we can compute the mean growth rate in the reactor

µ(s, I 0 (t) , x) = 1 L L 0 µ z (s, I (xz, t))dz,
where L is the depth of the reactor and where we have supposed that, even though the growth rate is not homogeneous in the reactor due to the light attenuation, the concentrations of s and x are kept homogeneous through continuous reactor stirring. It is this average growth rate that will be used in the lumped model that we develop.

We then have where µ is hypothetical maximum growth rate, a is a light attenuation coefficient, L is the depth of the reactor, K I represents the half-saturation coefficient relative to the light, K s represents the halfsaturation coefficient (Michaelis-Menten constant) relative to the substrate, ρ is the respiration rate and all the rest of parameters follows from (2.1). This model is shown (in a reduced notation form) in (2.2).

µ(s, I 0 (t) , x) = µ L L 0 I(xz,t) I(xz,t)+K I dz s s+K s = µ L L 0 I 0 (t)e

B Simplification and reduction of the system

It seems clear that the larger s translating into large growth rates. In a batch process, the initial substrate concentration should then always be very large so as to always keep the substrate in the region where s s+K s ≈ 1, i.e., at such levels that the growth rate only dependent of light influence and biomass at each instant.

Remark

The previous hypothesis, which allows a strong simplification of model, is used in [START_REF] Grognard | Optimization of a photobioreactor biomass production using natural light[END_REF] to change and reduce an system similar to (2.2) in order to consider in the dynamics of the biomass only the influence of light and thus reduce to a optimal control problem with an independent variable (biomass). In modelling terms, we assume that don't exist substrate limitation. This is always true in the initial times in a batch bioreactor, however, in long terms is not always true. By this reason high initial substrate concentrations are needed. We assume this hypothesis in our model. So, it is can be studied the reduced model

dx dt = µ axL ln I 0 (t)+K I I 0 (t)e -axL +K I x -ρx, ( 2.14)
which then encompasses all the relevant dynamics for the optimization problem. In order to more precisely determine the model, We must indicate what the varying light will be like. Classically, it is considered that day light varies as the square of a sinusoidal function so that

I 0 (t) = max sin 2πt T d , 0 2 ,
where T d is the length of the day. The introduction of such a varying light would however render the computations analytically untractable. Suppose that the light and dark periods appears in periodic intervals, for k = 0, 1, ..., and considering t 0 = 0 < t 1 < ... < t 2k < t 2k+1 < ..., we can approximate the light source by a step function:

I 0 (t) =    I 0 , t 2k ≤ t < t 2k+1 , ( light phase)
, 0, t 2k+1 ≤ t < t 2k+2 , (dark phase).

(2.15)

From (2.15) the biomass growth in the presence of light is reduced to

µ(I 0 (t) , x) =      µ axL ln I 0 +K I I 0 e -axL +K I , t 2k ≤ t < t 2k+1 , 0, t 2k+1 ≤ t < t 2k+2 .
Finally, we consider a last simplification to the model: instead of considering that the biomass growth in the presence of light has the form

µ 1 (x) = µ axL ln I 0 +K I I 0 e -axL +K I x,
which is an increasing and bounded function, following the idea in [START_REF] Grognard | Optimization of a photobioreactor biomass production using natural light[END_REF], we replace µ 1 (x) with another similar increasing bounded function given by

µ 2 (x) = µx c+x ,
where c is a fitting parameter (see appendix C).

From the simplifications above, the reduced system is

dx dt = µ(t)x c+x -ρx, ( 2.16)
where

µ (t) =    µ, t 2k ≤ t < t 2k+1 (light phase), 0, t 2k+1 ≤ t < t 2k+2 (dark phase).
(2.17)

C Parameter estimation

In order to use the model in a case study, we consider algae Chlamydomonas reinhardtii. The following table contains parameter values obtained in some references. For microalgae production is better to use a high incident light. This study was done under over-saturating photon flux density during the light period, I 0 = 700 [µmol. m -2 .s -1 ] [58]. Now, [START_REF] Kliphuis | Metabolic modeling of chlamydomonas reinhardtii: energy requirements for photoautotrophic growth and maintenance[END_REF] gives an respiration rate estimate of 0.53 [µmol O 2 . g -1 .h -1 ] and a molecular weight of 22.45 [g. molC -1 ] and from [START_REF] Falkowski | Aquatic photosynthesis[END_REF] the respiratory quotient of this mocroalgae is 1.15 [molC. molO To determine the light attenuation coefficient, we will use a simple linear regression with the data obtained from [START_REF] Kliphuis | Metabolic modeling of chlamydomonas reinhardtii: energy requirements for photoautotrophic growth and maintenance[END_REF]. Finally, it is necessary to estimate the last parameters by comparing the functions µ 1 (x) and µ 2 (x) that was defined above (appendix B) by a nonlinear least squares method. In the next table we shown the parameters were estimated by the nonlinear least squares method, from the data obtained in references (see table 2.5) for the model applied to C. reinhardtii with corresponding measurement units. Using the same procedure, we can obtain different parameters values for this model, for example by varying the incidental light I 0 (see table 2.3).

Chapter 3

Modelling and stability analysis of a microalgal pond with nitrification

Introduction

Microalgae culture is an emerging solution for a wide range of applications: food production, wastewater treatment, biofuel, among others. For large-scale production, microalgae can be fed with wastewater or digestate, that contains a large amount of ammonium [START_REF] García | High rate algal pond operating strategies for urban wastewater nitrogen removal[END_REF][START_REF] Uggetti | Anaerobic digestate as substrate for microalgae culture: the role of ammonium concentration on the microalgae productivity[END_REF]. The presence of oxygen (produced by photosynthesis) and ammonium favor the growth of nitrifiers, which transform ammonium into nitrate. On the other hand, microalgae can grow either on ammonium or nitrate as a nitrogen source. The conditions in which microalgae and nitrifiers can coexist is an important issue for microalgae pond operation.

The chemostat is a laboratory bioreactor in which fresh medium is continuously added and culture liquid is continuously removed, so that the culture volume remains constant. It is an ideal model for studying competition between exploited species. It is also used as a model of the wastewater treatment process. It has important industrial applications, including the commercial production of genetically modified organisms [START_REF] Smith | The Theory of the Chemostat[END_REF].

The basic mathematical model for the chemostat was first presented in [START_REF] Monod | La technique de la culture continue: Theorie et applications[END_REF][START_REF] Novick | Description of the chemostat[END_REF]. However, different proposals for mathematical modelling of the dynamics of multiple species in the chemostat have been made over the recent years, both for one limiting substrate (see [START_REF] Butler | Coexistence of competing predators in a chemostat[END_REF][START_REF] Hsu | A Mathematical Theory for Single-Nutrient Competition in Continuous Cultures of Micro-organims[END_REF][START_REF] Lobry | A new hypothesis to explain the coexistence of n species in the presence of a single resource[END_REF][START_REF] Mazenc | Stability analysis of an ecological model[END_REF][START_REF] Wolkowicz | Global Dynamics of a Mathematical Model of Competition in the Chemostat: General Response Functions and Differential Death Rates[END_REF]) and for more than one substrate (see [4,[START_REF] Ballyk | Exploitative Competition in the Chemostat for Two Perfectly Substitutable Resources[END_REF][START_REF] León | Competition between two species for two complementary or substitutable resources[END_REF][START_REF] Mazenc | Stability and stabilization for models of chemostats with multiple limiting substrates[END_REF]).

In the case of competition for one resource [START_REF] Grover | Resource Competition[END_REF][START_REF] Tilman | Resource competition and community structure[END_REF], for a well-mixed continuous culture (e.g. This chapter is based on the paper F. Mairet, H. Ramírez C. and A. Rojas-Palma, Modelling and stability analysis of a microalgal pond with nitrification, submitted to Applied Mathematical Modelling, Elsevier (2016). the chemostat), the competitive exclusion principle [5] provides conditions on growth rates under which only one species can generically persist. This has led to a large literature (for example, [5,[START_REF] Gajardo | Competitive exclusion principle in a model of chemostat with delays[END_REF][START_REF] Lobry | Effect on persistence of intra-specific competition in competition models[END_REF] and references therein) whose aim is to explain the discrepancy between the competitive exclusion principle and the fact that it is common for multiple competing species to survive in nature on one limiting substrate [START_REF] Mazenc | Stability and stabilization for models of chemostats with multiple limiting substrates[END_REF].

The case of two species and two substrates has been well studied in [START_REF] Ballyk | Global analysis of competition for perfectly substitutable resources with linear response[END_REF][START_REF] Ballyk | Exploitative Competition in the Chemostat for Two Perfectly Substitutable Resources[END_REF][START_REF] Macarthur | Species packing and competitive equilibrium for many species[END_REF]. In particular, these papers focus on ensuring the coexistence of two species based on the consideration of nullclines, the locations and stability of single species equilibria, and global stability results for coexistence are given. Nonetheless, the particular case of microalgae and nitrifiers -together with competition and cross-feeding -has never been addressed to the best of our knowledge.

In this work, we propose a reduced model of a microalgal pond with nitrification. This model considers the competition for nitrogen between microalgae and nitrifiers, but others possible interactions (CO2, O2) are neglected. Nitrification is represented by a one-step bioreaction. Microalgae can grow either on ammonium or nitrate, with a preference for ammonium. Light limitation (by self-shading) is also included in the microalgae growth rate.

The goal of this manuscript is to propose a dynamical model of the competition in a microalgal pond and to mathematically analyze this system. We finally provide some hints on how to manage the presence of nitrifiers.

Preliminaries.

1) By positive equilibrium point, we refer to an equilibrium point whose components are positive.

ii) Throughout this paper, global asymptotic stability should be interpreted as global asymptotic stability in the defined domain of the considered system, i.e., where the initial conditions have sense in terms of their meaning in the studied problem (For instance, they should be positives)

iii) By equilibrium existence, we mean that equilibrium points belong to the domain of definition of the considered system. iv) For convenience, we sometimes consider x ∈ R n , and for f : R n → R, we denote ∇ f = (

∂ f ∂x 1 , ..., ∂ f ∂x n ).

Problem statement

In first place, we define in Table 3.1 the state variables used in the algal pond model.

Stoichiometric Equations

We consider the following reactions:

1. Nitrifying Biomass Growth:

k 1 s 1 µ 1 (s 1 )x 1 -→ x 1 + k 1 s 2 . (3.1)
Note that we assume that all the nitrogen consumed is oxidized into nitrate, neglecting the nitrogen in the biomass. This assumption can be checked using ASM1 parameter values [1]: the nitrogen content of the biomass i XB (0.086 gN/gCOD) is negligible compared to consumption yield Y -1 A (4.17 gN/gCOD).

2. Microalgal Growth on Ammonium:

k 2 s 1 µ 2 (s 1 ,x 2 )x 2 -→ x 2 . (3.2)
3. Microalgal Growth on Nitrate:

k 2 s 2 µ 3 (s 1 ,s 2 ,x 2 )x 2 -→ x 2 . (3.3)

Kinetic Equations

Assumptions 3.1 We give some assumptions about growth functions of the different species involved in the dynamics.

H1

The function µ 1 (s 1 ) is iqual to zero for s 1 = 0 (onwards zero at zero), positive, continuous, increasing and bounded function (a Monod growth fulfills these assumptions for instance).

H2

The function µ 2 (s 1 , x 2 ) is defined by

µ 2 (s 1 , x 2 ) = μ2 (s 1 )φ 2 (x 2 ), (3.4)
where μ2 (s 1 ) is a zero at zero, positive, differentiable, increasing and bounded function and the function φ(x 2 ) is positive, differentiable, non increasing, φ 2 (0) = 1 and φ 2 (x 2 )x 2 is increasing and bounded (it can represent competitive inhibition).

H3

The function µ 3 (s 1 , s 2 , x 2 ) is defined by

µ 3 (s 1 , s 2 , x 2 ) = μ3 (s 2 )ϕ(s 1 )φ 2 (x 2 ), (3.5)
where μ3 (s 2 ) is a zero at zero, positive, differentiable, increasing and bounded function and the function ϕ(s 1 ) is positive, differentiable, non increasing, ϕ(0) = 1 and ϕ(s 1 )s 1 is increasing and bounded.

In particular we shall consider | μ3 (s 2 )| ≤ c|s 2 | for some c > 0 (as in Monod case) and in general

|µ 3 (x 2 , s 1 , s 2 )| ≤ c (x 2 , s 1 , s 2 ) .
The functions µ 1 (s 1 ), μ2 (s 1 ), and μ3 (s 2 ) are classical kinetics for substrate limitation. The function φ(x 2 ) is used to represent microalgae self-shading: the more biomass there is, the more light is attenuated and so the growth rate decreases. Finally, the function ϕ(s 1 ) is used to represent microalgae preference for ammonium over nitrate (i.e. ammonium inhibits growth on nitrate).

Mass Balance Equations

In a continuous reactor, the model is described as follow:

       ẋ1 = (µ 1 (s 1 ) -D)x 1 ẋ2 = (µ 2 (s 1 , x 2 ) + µ 3 (s 1 , s 2 , x 2 ) -D)x 2 ṡ1 = D(s in -s 1 ) -k 1 µ 1 (s 1 )x 1 -k 2 µ 2 (s 1 , x 2 )x 2 ṡ2 = -Ds 2 + k 1 µ 1 (s 1 )x 1 -k 2 µ 3 (s 1 , s 2 , x 2 )x 2 (3.6)
where D and s in are respectively the dilution rate and the ammonium input concentration. The system (3.6) is defined in the region

Ω = {(x 1 , x 2 , s 1 , s 2 ) ∈ R 4 | x 1 , x 2 , s 1 , s 2 ≥ 0}.
and all the parameters are positive. This model of competition in a chemostat is similar to the one presented in [START_REF] Ballyk | Global analysis of competition for perfectly substitutable resources with linear response[END_REF][START_REF] Ballyk | Exploitative Competition in the Chemostat for Two Perfectly Substitutable Resources[END_REF][START_REF] Hsu | Exploitative Competition of Microorganisms for Two Complementary Nutients in Continuous Cultures[END_REF] but considering intra-specific competition phenomenon through densitydependent growth functions [START_REF] Lobry | A new hypothesis to explain the coexistence of n species in the presence of a single resource[END_REF][START_REF] Lobry | Effect on persistence of intra-specific competition in competition models[END_REF][START_REF] Mazenc | Stability analysis of an ecological model[END_REF] and also cross-feeding (the nitrate produced by the nitrifiers can be consummed by the algae).

Model analysis

We consider in this section the mathematical analysis of the system (3.6). The first result is quite technical but important, because it implies that the system is well posed.

Lemma 3.2

For initial conditions in Ω, the solutions of the system (3.6) remain positive and are bounded for all t > 0.

PROOF. First of all, we need to prove that Ω is an invariant set. Suppose that s 1 = 0, where in this case ṡ1 = D(s in ) ≥ 0, then s 1 (t) ≥ 0, ∀t > 0. In the same way, if s 2 = 0 then ṡ2 = k 1 µ 1 (s 1 )x 1 ≥ 0, ∀t > 0. Finally, if x i = 0 for i = 1, 2, where ẋi = 0, all solution with initial condition x i (t 0 ) = 0 remains in this plane and the Picard-Lindelöf theorem implies that a solution in int(Ω) cannot cross these planes. So, each solution with initial condition in Ω remains in this set, i.e., this is an invariant set.

The dynamical behaviour of the solutions x 1 (t) of the first equation of (3.6) depends of s 1 (t) whose dynamics are described in the third equation of (3.6). We note (using comparision principle [START_REF] Khalil | Nonlinear Systems: Pearson New International Edition. Always learning[END_REF]) that

ṡ1 = D(s in -s 1 ) -k 1 µ 1 (s 1 )x 1 -k 2 µ 2 (s 1 , x 2 )x 2 ≤ D(s in -s 1 ) -k 1 µ 1 (s 1 )x 1 ,
then, we can consider the subsystem

ẋ1 = (µ 1 (s 1 ) -D)x 1 ṡ1 = D(s in -s 1 ) -k 1 µ 1 (s 1 )x 1 (3.7)
and defining r = k 1 x 1 + s 1 then we have

ṙ = k 1 ẋ1 + ṡ1 = D (s in -k 1 x 1 -s 1 ) = D (s in -r) ,
the following differential equation is obtained ṙ + Dr = Ds in , where r(t) = s in -(s in -r(0))e -Dt ≤ s in , ∀r(0) ≤ s in . Otherwise, if r(0) > s in then r(t) = s in + (r(0) -s in )e -Dt is decreasing and its maximum value is r(0). If we denote L = max{r(0), s in } then r(t) ≤ L, ∀t > 0. In particular, k 1 x 1 ≤ L in this case.

To show that all solutions are bounded, we consider w = k 1 x 1 + k 2 x 2 + s 1 + s 2 then (using majorization)

ẇ = k 1 ẋ1 + k 2 ẋ2 + ṡ1 + ṡ2 = D (s in -k 1 x 1 -k 2 x 2 -s 1 -s 2 ) + µ 1 (s 1 )k 1 x 1 = D (s in -w) + µ 1 (s 1 )k 1 x 1 ≤ D (s in -w) + µ 1 (s in )L ≤ Ds in + µ 1 (s in )L -Dw, the following inequality is obtained 0 ≤ ẇ + Dw ≤ Ds in + µ 1 (s in )L.
From a theorem on differential inequalities [START_REF] Szarski | Differential Inequalities[END_REF]

0 ≤ w(t) ≤ s in + L D µ 1 (s in ) + (w(0) -(s in + L D µ 1 (s in )
))e -Dt , ∀t > 0 when t → ∞ we get w ≤ s in + L D µ 1 (s in ). Therefore, by defining the set

R = (x 1 , x 2 , s 1 , s 2 ) ∈ Ω | k 1 x 1 + k 2 x 2 + s 1 + s 2 < max{w(0), s in + L D µ 1 (s in )} ,
the set R is the region where all solution trajectories of the system (3.6) with initial conditions in Ω is confined.

Remark

1. The dynamical system is said to be dissipative if all positive trajectories eventually lie in a bounded set [START_REF] Smith | The Theory of the Chemostat[END_REF]. This is sufficient to ensure that all solutions exist for all positive time. The last lemma ensure this property for the system (3.6).

2. In the first equation in (3.6) using the comparision principle [START_REF] Khalil | Nonlinear Systems: Pearson New International Edition. Always learning[END_REF]

ẋ1 = (µ 1 (s 1 ) -D) x 1 ≤ (µ 1 (s in ) -D) x 1 , we obtain x 1 (t) ≤ x 1 (0)e -(D-µ 1 (s in ))t ,
Suppose that the condition µ 1 (s in ) < D is fulfilled. Then x 1 → 0 when t → ∞ ∀x 1 (0) > 0, i.e., if the dilution rate is greater than the growth rate of nitrifying biomass at input ammonium level, their concentration tends to disappear in time. In the same way, we can prove a similar result for x 2 (t) considering in the second equation of (3.6) that s 2 ≤ s in + L D µ 1 (s in ) = s2 and µ 2 (s in , 0) + µ 3 (s in , s2 , 0) < D, i.e., the biomass concentration of microalgae tends to zero in this case.

Equilibrium existence and local stability Lemma 3.3 (non-coexistence equilibrium points)

The non coexistence equilibrium points of the system (3.6) are 1. (Washout) E 0 = (0, 0, s in , 0) which always exists.

(Nitrifier only) E

n = ( s in -λ n k 1 , 0, λ n , s in -λ n ) with λ n = µ -1 1 (D)
which exists if and only if D < µ 1 (s in ) (or equivalently λ n < s in ).

(Microalgae only

) E a = (0, s in -λ s k 2 , λ s , 0) with λ s the unique solution of µ 2 (s 1 , s in -s 1 k 2 ) = D which exists if and only if D < µ 2 (s in , 0) (or equivalently λ s < s in ).

PROOF. The equilibrium points are the positive solutions of the system of nonlinear equations

(µ 1 (s 1 ) -D)x 1 = 0, (µ 2 (s 1 , x 2 ) + µ 3 (s 1 , s 2 , x 2 ) -D)x 2 = 0, D(s in -s 1 ) -k 1 µ 1 (s 1 )x 1 -k 2 µ 2 (s 1 , x 2 )x 2 = 0, -Ds 2 + k 1 µ 1 (s 1 )x 1 -k 2 µ 3 (s 1 , s 2 , x 2 )x 2 = 0. (3.8)
From the first two equations, we have In this setting μn = 2, μs = 1, K n = 1, K s = 0.12, K x = 0.5, s in = 1, D = 0.5, k 1 = 1 and considering the growth functions in the application section, the graph show the existence of the three equilibrium points that was proven in lemma above. Now, if any other equilibrium point of (3.6) exists, there must be a positive equilibrium. Assuming this fact, the positive solutions of equation (3.8) must simultaneously solve the equations

µ 1 (s 1 ) = D, μ2 (s 1 ) + μ3 (s 2 )ϕ(s 1 ) = x 2 φ(x 2 ) D, (3.9)
D(s in -s 1 ) = k 1 µ 1 (s 1 )x 1 + k 2 μ2 (s 1 )v(x 2 ), -Ds 2 = -k 1 µ 1 (s 1 )x 1 + k 2 μ3 (s 2 )ϕ(s 1 )v(x 2 ), (3.10)
where v(x 2 ) = x 2 φ(x 2 ). If we denote y = v(x 2 ), using Cramer's rule, a solution (x * 1 , x * 2 ) = (x * 1 , v -1 (y * )) of system (3.10) for s 1 , s 2 fixed is given by

x * 1 = k 2 D((s in -s 1 ) μ3 (s 2 )ϕ(s 1 )-s 2 μ2 (s 1 )) d(s 1 ,s 2 ) , v(x * 2 ) = k 1 Dµ 1 (s 1 )(s in -s 1 -s 2 ) d(s 1 ,s 2 )
, where d(s 1 , s 2 ) = k 1 k 2 µ 1 (s 1 )( μ2 (s 1 ) + μ3 (s 2 )ϕ(s 1 )) > 0, ∀s 1 , s 2 > 0.

When s 1 and s 2 are positive solutions of (3.9) then we obtain the result of proposition 3.4. In the following, we define

f (s 1 , x 2 ) = µ 2 (s 1 , x 2 ) + µ 3 (s 1 , s in -s 1 -k 2 x 2 , x 2 ).
(3.11)

Proposition 3.4 (coexistence)

A necessary and sufficient condition for the existence and uniqueness of a positive equilibrium point of (3.6) is given by D < min {µ 1 (s in ), µ 2 (s in , 0), µ 2 (λ n , 0) + µ 3 (λ n , s inλ n , 0)} (3.12)

or equivalently λ n < s in , λ s < s in , and D < γ n,0 , with

γ n,0 := f (λ n , 0) = µ 2 (λ n , 0) + µ 3 (λ n , s in -λ n , 0).
Furthermore, in this case, the unique positive equilibrium E

c = (x * 1 , x * 2 , s * 1 , s * 2 ) is described by s * 1 = λ n , s * 2 = s in -λ n -k 2 x * 2 , x * 1 = s in -λ n -k 2 x * 2 k 1 + k 2 µ 3 (λ n ,x * 2 )x * 2 Dk 1 ,
where x * 2 is a positive solution of the nonlinear equation

µ 2 (λ n , x 2 ) + µ 3 (λ n , s in -λ n -k 2 x 2 , x 2 ) = D.
PROOF. Suppose that µ 1 (s 1 ) = D and µ 2 (s 1 , If we add the last two equations in (3.8) then

x 2 ) + µ 3 (s 1 , s 2 , x 2 ) = D.
D(s in -s 1 -s 2 ) -k 2 (µ 2 (s 1 , x 1 ) + µ 3 (s 1 , s 2 , x 2 ))x 2 = 0 D(s in -s 1 -s 2 -k 2 x 2 ) = 0,
that is, the equilibrium points must necessarily fulfill the equality relation

s in -s 1 -s 2 -k 2 x 2 = 0. (3.13) then necessarily s 2 = s in -λ n -k 2 x 2 .
From the third equation of (3.8), using the equality

µ 1 (s 1 ) = D D(s in -λ n ) -k 1 Dx 1 -k 2 µ 2 (λ n , x 2 )x 2 = 0,
therefore, solving for the variable x 1 we obtain

x * 1 = D(sin-λ n )-k 2 µ 2 (λ n ,x 2 )x 2 Dk 1 , or equivalently x * 1 = s in -λ n -k 2 x 2 k 1 + k 2 µ 3 (λ n ,x 2 )x 2 Dk 1 .
Finally, x * 2 should be the solution of the equation

µ 2 (λ n , x 2 ) + µ 3 (λ n , s * 2 , x 2 ) = D (3.14)
It is noted that by replacing the term s * 2 into the last equation, the unique variable is x 2 . From assumptions H1-H3 the function (3.11) is positive in ]0, ∞[×]0, ∞[, increasing with respect to s 1 and is strictly decreasing with respect to x 2 . We note that f (s in , 0) = µ 2 (s in , 0) so a necessary condition for the existence of the positive equilibrium is D < µ 2 (s in , 0). Finally, the maximum possible value of the function evaluated in the equilibrium value λ n is f (λ n , 0) = µ 2 (λ n , 0) + µ 3 (λ n , s inλ n , 0) and then to guarantee the existence of solution of the equation (3.14)

f (λ n , 0) = µ 2 (λ n , 0) + µ 3 (λ n , s in -λ n , 0) := γ n,0 > D.
this completes the proof.

Table 3.2. Equilibrium existence conditions for break-even concentrations

Equilibria Existence

E 0 Always E n λ n < s in E a λ s < s in E c λ n < s in , λ s < s in , D < γ n,0
Proposition 3.5 (Local stability of non-coexistence equilibrium points)

1. The equilibrium point E 0 = (0, 0, s in , 0) is locally asymptotically stable if µ 1 (s in ) < D and µ 2 (s in , 0) < D (or equivalently s in < λ n and s in < λ s ).

The equilibrium point E

n = ( s in -λ n k 1 , 0, λ n , s in -λ n ) is locally asymptotically stable if µ 2 (λ n , 0) + µ 3 (λ n , s in -λ n , 0) < D (or equivalently γ n,0 < D).
3. The equilibrium point E a = (0, s in -λ s k 2 , λ s , 0) is locally asymptotically stable if µ 1 (λ s ) < D (or equivalently λ s < λ n ).

PROOF. We shall use the Lyapunov indirect method [START_REF] Chicone | Ordinary Differential Equations with Applications[END_REF].

1. The Jacobian matrix of system (3.6) (see (3.24) in appendix A) evaluated in the point E 0

J(E 0 ) =     µ 1 (s in ) -D 0 0 0 0 µ 2 (s in , 0) -D 0 0 -k 1 µ 1 (s in ) -k 2 µ 2 (s in , 0) -D 0 k 1 µ 1 (s in ) 0 0 -D    
is a lower triangular matrix. Then its eigenvalues are the diagonal elements. In this case, all the eigenvalues are negatives (negative real part) if the conditions µ 1 (s in ) < D and µ 2 (s in , 0) < D are satisfied.

2. The Jacobian matrix of the system (3.6) evaluated in the point E n is

J(E n ) =      0 0 s in -λ n k 1 dµ 1 ds 1 (λ n ) 0 0 J 22 0 0 -k 1 D -k 2 µ 2 (λ n , 0) J 33 0 k 1 D -k 2 µ 3 (λ n , s in -λ n , 0) (s in -λ n ) dµ 1 ds 1 (λ n ) -D     
where µ 1 is increasing,

dµ 1 ds 1 (λ n ) > 0 and J 22 = µ 2 (λ n , 0) + µ 3 (λ n , s in -λ n , 0) -D = 0, J 33 = -D -(s in -λ n ) dµ 1 ds 1 (λ n ) < 0.
The eigenvalues of the matrix J(E n ) are {-D, J 22 } and the roots of the quadratic equation

λ 2 -J 33 λ + (s in -λ n ) dµ 1 ds 1 (λ n )D = 0, whose discriminant is △ = (J 33 ) 2 -4(s in -λ n ) dµ 1 ds 1 (λ n )D = ((s in -λ n ) dµ 1 ds 1 (λ n ) -D) 2 > 0,
i.e., the roots of the last equation are real and (from Descartes' rules of signs) both roots are negative, then if µ 2 (λ n , 0) + µ 3 (λ n , s inλ n , 0) < D the equilibrium point E n is locally asymptotically stable.

3. The Jacobian matrix of the system (3.6) evaluated in the point E a is

J(E a ) =     µ 1 (λ s ) -D 0 0 0 0 J 22 J 23 s in -λ s k 2 ∂µ 3 ∂s 2 (λ s , 0, s in -λ s k 2 ) -k 1 µ 1 (λ s ) J 32 J 33 0 k 1 µ 1 (λ s ) 0 0 -D    
where

J 22 = s in -λ s k 2 ∂µ 2 ∂x 2 (λ s , s in -λ s k 2 ) < 0, J 23 = s in -λ s k 2 ∂µ 2 ∂s 1 (λ s , s in -λ s k 2 ) > 0, J 32 = -(s in -λ s ) ∂µ 2 ∂x 2 (λ s , s in -λ s k 2 ) -k 2 D = 0, J 33 = -D -(s in -λ s ) ∂µ 2 ∂s 1 (λ s , s in -λ s k 2 ) < 0.
The eigenvalues of the matrix J(E a ) are {µ 1 (λ s ) -D, J 22 } (where µ 1 (λ s ) -D < 0 when λ s < λ n ) and the eigenvalues of the submatrix

A = J 22 J 23 J 32 J 33 ∈ M 2×2 (R),
determined by the equation

λ 2 -trace(A)λ + det(A) = 0, (3.15)
where trace(A) = J 22 + J 33 and det(A) = J 22 J 33 -J 23 J 32 . In this case

trace(A) = s in -λ s k 2 ( ∂µ 2 ∂x 2 (λ s , s in -λ s k 2 ) -k 2 ∂µ 2 ∂s 1 (λ s , s in -λ s k 2 )) -D, det(A) = -s in -λ s k 2 ( ∂µ 2 ∂x 2 (λ s , s in -λ s k 2 ) -k 2 ∂µ 2 ∂s 1 (λ s , s in -λ s k 2 ))D.
The discriminant of (3.15) is

△ = (trace(A)) 2 -4 det(A) = ( s in -λ s k 2 ( ∂µ 2 ∂x 2 (λ s , s in -λ s k 2 ) -k 2 ∂µ 2 ∂s 1 (λ s , s in -λ s k 2 )) + D) 2 > 0,
then the roots of (3.15) are real and are given by

λ 1 = -2D < 0 λ 2 = 2 s in -λ s k 2 ( ∂µ 2 ∂x 2 (λ s , s in -λ s k 2 ) -k 2 ∂µ 2 ∂s 1 (λ s , s in -λ s k 2 )) < 0,
and where all the eigenvalues of the matrix J(E a ) are negative (negative real part) and none is equal to zero. The equilibrium point E a is locally asymptotically stable when λ s < λ n .

From the three cases above, the statement is proved.

E 0 s in < λ n , s in < λ s E n γ n,0 < D E a λ s < λ n
In the next lemma we investigate the nature of the coexistence equilibrium point.

Lemma 3.6 Suppose the conditon (3.12), the equilibrium E c is an hyperbolic equilibrium point.

PROOF. The determinant of Jacobian matrix of the system (3.6) (see (3.24)) following the same notation from the previous proposition, is given by

det(J(E c )) = k 1 x * 1 D dµ 1 ds 1 (λ n ) (J 22 J 44 -J 32 J 24 -J 24 J 42 ) = x * 2 D k 2 ∂µ 3 ∂s 2 (λ n , s * 2 , x * 2 ) -∂µ 2 ∂x 2 (λ n , x * 2 ) + ∂µ 3 ∂x 2 (λ n , s * 2 , x * 2 ) .
Suppose parameters are varied so that an eigenvalue of E c passes through zero. At the bifurcation point, either E c collapses with the equilibrium points in ∂Ω (the boundary of the set Ω) or it exists as an equilibrium in int(Ω), in which case x * i > 0 and s * i > 0 (i = 1, 2) then the determinant of E c cannot be zero (it is positive, in fact it coincides with the coefficient a 4 of the characteristic polynomial as shown in the previous proposition) and so there is no bifurcation. Thus, if parameters are varied so that an eigenvalue of J(E c ) reaches zero, then E c collapses either with E n or E a (or possibly E 0 if there is a higher order bifurcation). The next proposition is a direct consequence of the results in Section 3.1 for the reduced system (3.18) Proposition 3.7 Let us consider the system (3.18) a) For initial conditions in Ω 2 , the solutions of the system remain positive and are uniformly bounded for all t > 0 in the region

R 2 = (x 1 , x 2 , s 1 ) ∈ Ω 2 | k 1 x 1 + k 2 x 2 + s 1 < max{w(0), L D µ 1 (s in )} .
b) The equilibrium points of the system are the following ones 1. Ē0 = (0, 0, s in ) that always exists.

Ēn

= ( s in -λ n k 1 , 0, λ n ) with λ n = µ -1 1 (D) that exists if and only if D < µ 1 (s in ). 3. Ēa = (0, s in -λ s k 2 , λ s ) with λ s the unique solution of µ 2 (s 1 , s in -s 1 k 2 ) = D that exists if and only if D < µ 2 (s in , 0). 4. Ēc = s in -λ n -k 2 x * 2 k 1 + k 2 µ 3 (λ n ,x * 2 )x * 2 Dk 1 , x * 2 , λ n where x * 2 is a positive solution of equation µ 2 (λ n , x 2 ) + µ 3 (λ n , x 2 ) = D,
that exists and is unique if and only if

D < min {µ 1 (s in ), µ 2 (s in , 0), µ 2 (λ n , 0) + µ 3 (λ n , 0)} .
c) Local stability of non-coexistence equilibrium points 1. If µ 1 (s in ) < D and µ 2 (s in , 0) < D, then the equilibrium point Ē0 is locally asymptotically stable.

2. If µ 2 (λ n , 0) + µ 3 (λ n , 0) < D, then the equilibrium point Ēn is locally asymptotically stable.

3. If µ 1 (λ s ) < D, then the equilibrium point Ēa is locally asymptotically stable.

We are now able to show a sufficient condition for the local stability of the interior positive equilibrium point Ēc Proposition 3.8 Asuming condition (3.12), if

0 < δγ n,x * 2 , (3.19)
the equilibrium point Ēc is locally asymptotically stable, where

δγ n,x * 2 = ∂µ 2 ∂s 1 (λ n , x * 2 ) + ∂µ 3 ∂s 1 (λ n , x * 2 ). (3.20)
PROOF. The Jacobian matrix evaluated in the positive equilibrium point Ēc is

J( Ēc ) =   0 0 J 13 0 J 22 J 23 -k 1 D J 32 J 33   ,
where

J 13 = x * 1 dµ 1 ds 1 (λ n ) > 0, J 22 = ( ∂µ 2 ∂x 2 (λ n , x * 2 ) + ∂µ 3 ∂x 2 (λ n , x * 2 ))x * 2 < 0, J 23 = ( ∂µ 2 ∂s 1 (λ n , x * 2 ) + ∂µ 3 ∂s 1 (λ n , x * 2 ))x * 2 = 0, J 32 = -k 2 (x * 2 ∂µ 2 ∂x 2 (λ n , x * 2 ) + µ 2 (λ n , x * 2 )) < 0, J 33 = -D -k 1 x * 1 dµ 1 ds 1 (λ n ) -k 2 x * 2 ∂µ 2 ∂s 1 (λ n , x * 2 ) < 0,
The characteristic polynomial asociated to J( Ēc ) is

p(λ) = λ 3 + a 1 λ 2 + a 2 λ + a 3 ,
where

a 1 = -J 22 -J 33 > 0, a 2 = J 22 J 33 -J 13 J 31 -J 23 J 32 , a 3 = -k 1 DJ 13 J 22 > 0.
The Routh-Hurwitz criterion in this case provides that the roots of the polynomial p(λ) have negative real part if and only if a i > 0 (i = 1, 2, 3) and a 1 a 2 -a 3 > 0. Suppose that 0 <

∂µ 2 ∂s 1 (λ n , x * 2 ) + ∂µ 3 ∂s 1 (λ n , x * 2 )
we obtained that J 23 > 0, with this a 2 > 0 and

a 1 a 2 -a 3 = J 13 J 31 J 33 -J 2 22 J 33 + J 22 J 2 33 + J 23 J 32 (J 22 + J 33 ) > 0,
then the eigenvalues of J( Ēc ) have negative real part, i.e., the coexistence equilibrium point Ēc is locally asymptotically stable.

s in < λ n , s in < λ s Ēn λ n < s in γ n,0 < D Ēa λ s < s in λ s < λ n Ēc λ n < s in , 0 < δγ n,x * 2 λ s < s in , D < γ n,0

About the global behaviour

In Proposition 3.8 we prove that the coexistence equilibrium point of the system (3.18) can be locally asymptotically stable. To try to demonstrate the global stability of this point, a first approach would be to use the direct method of Lyapunov [START_REF] Chicone | Ordinary Differential Equations with Applications[END_REF][START_REF] Khalil | Nonlinear Systems: Pearson New International Edition. Always learning[END_REF], which gives sufficient conditions for global asymptotic stability. However, it is not always possible to construct such functions. For this reason, some researchers have used different approaches, for example in [START_REF] Ballyk | Exploitative Competition in the Chemostat for Two Perfectly Substitutable Resources[END_REF][START_REF] Vargas-De-León | Global stability of a population dynamics model with inhibition and negative feedback[END_REF] an non-strict Lyapunov function [START_REF] Mazenc | Stability and stabilization for models of chemostats with multiple limiting substrates[END_REF] is constructed and LaSalle invariance principle [START_REF] Sastry | Nonlinear Systems: Analysis, Stability, and Control[END_REF] is applied, in fact [START_REF] Ballyk | Exploitative Competition in the Chemostat for Two Perfectly Substitutable Resources[END_REF] give a different focus, using second additive compound matrices [START_REF] Li | On r.a. smith's autonomous convergence theorem[END_REF][START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF] for linear growth function in the case of exploitative competition for two resources. In [START_REF] Mazenc | Stability analysis of an ecological model[END_REF] for a model of n species in competition for a single resource, the authors make a change of coordinates and they construct a Lyapunov function that is not radially unbounded. However, they study the sign properties of the time derivative of the function and Barbalat's lemma [START_REF] Khalil | Nonlinear Systems: Pearson New International Edition. Always learning[END_REF] is applied to conclude the attractiveness of the positive equilibrium. For general chemostat models, in [START_REF] Hsu | A survey of constructing Lyapunov functions for mathematical models in population biology[END_REF] the author gives a complete survey about the construction of Lyapunov functions. However, in the same article some open problems are proposed. On other hand, in [START_REF] De Leenheer | Global stability in a chemostat with multiple nutrients[END_REF][START_REF] Hsu | A survey of mathematical models of competition with an inhibitor[END_REF] the authors reduce their model and they study the system of lower dimension using the theory of asymptotically autonomous systems [START_REF] Smith | The Theory of the Chemostat[END_REF][START_REF] Strauss | On asymptotically autonomous differential equations[END_REF]. We shall develop the same idea in this section.

Remark Some important considerations are now emphasized i) The system (3.18) can be considered as a perturbation of the system

   ẋ1 = (µ 1 (s 1 ) -D)x 1 ẋ2 = (µ 2 (s 1 , x 2 ) -D)x 2 ṡ1 = D(s in -s 1 ) -k 1 µ 1 (s 1 )x 1 -k 2 µ 2 (s 1 , x 2 )x 2
(3.21)

A similar system, for n species but replacing µ 1 (s 1 ) by µ 1 (s 1 )θ(x 1 ) where θ(x 1 ) is strictly decreasing and θ(0) = 1, was studied in [START_REF] Lobry | A new hypothesis to explain the coexistence of n species in the presence of a single resource[END_REF][START_REF] Lobry | Effect on persistence of intra-specific competition in competition models[END_REF][START_REF] Lobry | Persistence in ecological models of competition for a single resource[END_REF]. Under certain assumptions it is possible to prove the existence and uniqueness of a positive equilibrium globally asymptotically stable. Although the system (3.21) does not fulfill one of these hypotheses (in this case θ(x 1 ) ≡ 1, then it is not strictly decreasing), [START_REF] Lobry | Effect on persistence of intra-specific competition in competition models[END_REF] provides a demonstration based on the Lyapunov direct method that includes the case of system (3.21).

ii) We note also that the stability results for the system (3.18) are valid for the system (3.21) (considering µ 3 (λ n , x * 2 ) = 0, i.e., vanishing perturbation [START_REF] Khalil | Nonlinear Systems: Pearson New International Edition. Always learning[END_REF]). In particular, the Jacobian matrix of system (3.21) evaluated on their unique positive equilibrium point

E * = ( s in -λ n k 1 - k 2 x * 2 k 1 , x * 2 , λ n ) with x * 2 solution of µ 2 (λ n , x 2 ) = D, is Hurwitz under the condition (3.19) in Proposition 3.8.
Proposition 3.9 (Some results of global stability for non-coexistence points) 1. Suppose that µ 1 (s in ) < D and µ 2 (s in , 0) < D, the equilibrium point Ē0 is globally asymptotically stable for all initial condition.

2. Suppose that µ 2 (λ n , 0) + µ 3 (λ n , 0) < D, the single-species equilibrium Ēn is globally asymptotically stable with respect to all solutions for which x 1 (t 0 ) > 0, and x 2 (t 0 ) = 0.

3. Suppose that µ 1 (λ s ) < D, the single-species equilibrium Ēa is globally asymptotically stable with respect to all solutions for which x 2 (t 0 ) > 0, and x 1 (t 0 ) = 0.

PROOF. 1. If the necessary conditions for the local stability of the washout equilibrium point Ē0 are fulfilled then (from Proposition 3.7) this is the unique equilibrium point in Ω 2 and the solutions are bounded in this case, i.e., the washout equilibrium point is globally asymptotically stable on this set.

2. We note that x 2 = 0 is invariant, then the dynamics for (x 1 (t), s 1 (t)) is determined by the reduced system ẋ1 = (µ 1 (s 1 ) -D)x 1 ṡ1 = D(s in -s 1 ) -k 1 µ 1 (s 1 )x 1 the last system under assumption H1 represent the classical general chemostat model for one species, where for any x 1 (t 0 ) > 0, we have λ n < s in implies x 1 (t) → s in -λ n k 1 (see theorem 3.2 in [START_REF] Smith | The Theory of the Chemostat[END_REF]).

3. We note that x 1 = 0 is invariant, then the dynamics for (x 2 (t), s 1 (t)) is determined by the reduced system ẋ2 = (µ 2 (s 1 ,

x 2 ) + µ 3 (s 1 , x 2 ) -D)x 2 ṡ1 = D(s in -s 1 ) -k 2 µ 2 (s 1 )x 1 (3.22)
the system (3.22) can be considered as a perturbation of system

ẋ2 = (µ 2 (s 1 , x 2 ) -D)x 2 ṡ1 = D(s in -s 1 ) -k 2 µ 2 (s 1 )x 1 (3.23)
The system (3.23) was studied in [START_REF] Lobry | Effect on persistence of intra-specific competition in competition models[END_REF] in particular, the authors show the global asymptotical stability of the unique positive equilibrium through the direct method of Lyapunov, this equilibrium matches with Ēa in this case and where µ 3 (λ s , x * ) = 0, ie in the equilibrium Ēa the perturbation vanish, the Lemma 9.1 in [START_REF] Khalil | Nonlinear Systems: Pearson New International Edition. Always learning[END_REF] assure that the equilibrium is globally asymptotically stable.

Theorem 3.10 (Ultimate bounded solutions near positive equilibrium) Assuming condition (3.19), the solutions of (3.18) with initial condition in R 2 tend asymptotically to the positive coexistence equilibrium point Ēc in a finite period of time and then remain bounded and near equilibrium.

PROOF. The demonstration will be given in two steps. First of all, we shall prove that the system (3.21) have a positive equilibrium point globally asymptotically stable trought the Lyapunov direct method, following the construction in [START_REF] Lobry | Effect on persistence of intra-specific competition in competition models[END_REF]. Secondly, we use a result of non-vanishing perturbated systems from [START_REF] Khalil | Nonlinear Systems: Pearson New International Edition. Always learning[END_REF] and the last remark to conclude the statement. The system (3.21) can be written as

   ẋ1 = ( μ1 (s 1 )φ 1 (x 1 ) -d 1 )x 1 ẋ2 = ( μ2 (s 1 )φ 2 (x 2 ) -d 2 )x 2 ṡ1 = f (s 1 ) -k 1 μ1 (s 1 )φ 1 (x 1 )x 1 -k 2 μ2 (s 1 )φ 2 (x 2 )x 2
where f (s 1 ) = D(s in -s 1 ), μ1 (.) = µ 1 (.), φ 1 (x 1 ) ≡ 1 and d 1 = d 2 = D. Clearly f (0) > 0 and from hypothesis H1-H2 φ i (x i ) are positive, decreasing, φ i (0) = 1 and

φ i (x i )x i is increasing for i = 1, 2.
There is a unique positive equilibrium point (x

* 1 , x * 2 , s * 1 ) ∈]0, ∞[ 2 ×]0, s in [ such that μi (s * 1 )φ i x * i = D and ∑ 2 i=1 k i μi (s * i )φ i (x * i )x * i = f (s * 1 )
, the functions μi are bounded, zero at zero, increasing and μ′ i (0) > 0. There is a positive function Λ and positive constants c i with i = 1, 2 such that

Λ(0) = c i μi (s * 1 ) μ′ i (s * 1 )s * 1 ,
and, for all s 1 > 0,

s 1 = s * 1 , Λ(s 1 ) = c i s 1 μi (s 1 )
μi (s 1 ) -μi (s * 1 )

s 1 -s * 1 .
Also, we can define a function

Γ(s 1 ) = - f (s 1 ) -f (s * 1 ) + ∑ 2 i=1 k i ( μi (s 1 ) -μi (s * 1 ))φ i (x * i )x * i s 1 -s * 1 ,
which is positive under the above assumptions. From [START_REF] Lobry | Effect on persistence of intra-specific competition in competition models[END_REF] we have that

V(x 1 , x 2 , s 1 ) = c 1 x 1 -x * 1 0 l l+x * 1 dl + c 2 x 2 -x * 2 0 l l+x * 2 dl +c 2 x * 2 x 2 -x * 2 0 φ 2 (l+x * 2 )-φ 2 (x * 2 ) φ 2 (l+x * 2 )
dl +

s 1 -s * 1 0 Λ(l+s * 1 )l l+s * 1 dl,
is a Lyapunov function for the system (3.21) which is positive definite and proper on

D t =] - x * 1 , ∞[×] -x * 2 , ∞[×] -s * 1 ,
∞[ (for details, see [START_REF] Lobry | Effect on persistence of intra-specific competition in competition models[END_REF]), and its derivative along the trajectories satisfies

V = -W(x 1 , x 2 , s 1 ), with W(x 1 , x 2 , s 1 ) = Λ(s 1 )Γ(s 1 ) s 1 (s 1 -s * 1 ) 2 + c 2 α(x 2 )β(x 2 ) φ 2 (x 2 ) (x 2 -x * 2 ) 2 ,
where α(x 2 ) = -µ 2 (s * 1 )

φ 2 (x 2 )-φ 2 (x * 2 ) x 2 -x * 2 and β(x 2 ) = φ 2 (x 2 )x 2 -φ 2 (x * 2 )x * 2 x 2 -x * 2
are positive from the assumptions. Then the positive equilibrium (x * 1 , x * 2 , s * 1 ) is a globally asymptotically and a locally exponentially stable equilibrium point of the system (3.21). Now we note that the Jacobian matrix of system (3.21) evaluated in their unique positive equilibrium point E * is Hurwitz under the condition (3.19) of Proposition 3.8, i.e., all their eigenvalues are strictly in the left half plane. Besides, if we denote

f 1 (x 1 , x 2 , s 1 ) =   ( μ1 (s 1 )φ 1 (x 1 ) -d 1 )x 1 ( μ2 (s 1 )φ 2 (x 2 ) -d 2 )x 2 f (s 1 ) -k 1 μ1 (s 1 )φ 1 (x 1 )x 1 -k 2 μ2 (s 1 )φ 2 (x 2 )x 2   ,
the vector field associated to system (3.21) and

f 2 (x 1 , x 2 , s 1 ) =   0 µ 3 (s 1 , x 2 )x 2 0   ,
the perturbation vector, and considering

∇V(x 1 , x 1 , s 1 ) = c 1 x 1 -x * 1 x 1 , c 2 x 2 -x * 2 x 2 -c 2 x * 2 φ 2 (x 2 )-φ 2 (x * 2 ) φ 2 (x 2 )
, Λ(s 1 )

s 1 -s * 1 s 1 ,
it is straightforward (see Lemma 4.3 in [60]) to note that there exist constants C i with i = 1, 2, 3 in terms of the parameters of the system, such that 2 , and from hypothesis H3, the non-vanishing perturbation term is increasing and bounded

C 1 (x 1 -x * 1 , x 2 -x * 2 , s 1 -s * 1 ) ≤ V(x 1 , x 2 , s 1 ) ≤ C 2 (x 1 -x * 1 , x 2 -x * 2 , s 1 -s * 1 ) , V = ∂V ∂t + ∇V • f 1 (x 1 , x 2 , s 1 ) = -W(x 1 , x 2 , s 1 ) ≤ -C 3 (x 1 -x * 1 , x 2 -x * 2 , s 1 -s * 1 )
|µ 3 (s 1 , x 2 )| = | μ3 (s in -k 2 x 2 -s 1 )||ϕ(s 1 )||φ 2 (x 2 )| ≤ | μ3 (s in -k 2 x 2 -s 1 )| ≤ μ3 (s in ). Now ∇V • f 2 (x 1 , x 2 , s 1 ) = c 2 x 2 -x * 2 x 2 -c 2 x * 2 φ 2 (x 2 )-φ 2 (x * 2 ) φ 2 (x 2 ) µ 3 (s 1 , x 2 )x 2 = c 2 (x 2 -x * 2 ) -x * 2 x 2 φ 2 (x 2 )-x 2 φ 2 (x * 2 ) φ 2 (x 2 ) µ 3 (s 1 , x 2 ),
from hypothesis H2-H3 and using the mean value theorem, there exists ξ ∈ [x 2 , x * 2 ] such that

x 2 φ 2 (x 2 ) -x 2 φ 2 (x * 2 ) = x 2 φ 2 (x 2 ) -x * 2 φ 2 (x * 2 ) + x * 2 φ 2 (x * 2 ) -x 2 φ 2 (x * 2 ) = (ξφ 2 (ξ)) ′ (x 2 -x * 2 ) -φ 2 (x * 2 )(x 2 -x * 2 ) = (ξφ 2 (ξ)) ′ -φ 2 (x * 2 ))(x 2 -x * 2 ),
by using this equality in the above expression, we have

∇V • f 2 (x 1 , x 2 , s 1 ) = c 2 φ 2 (x 2 ) -x * 2 ((ξφ 2 (ξ)) ′ -φ 2 (x * 2 )) μ3 (s in -k 2 x 2 -s 1 )ϕ(s 1 )(x 2 -x * 2 ),

Application

In [1], a dynamical model of an algal pond for wastewater treatment is proposed. In line with this model, let us consider the Monod growth functions

µ 1 (s 1 ) = μn s 1 K n + s 1 , μ2 (s 1 ) = μs s 1 K s + s 1 , μ3 (s 2 ) = μp s 2 K p + s 2 ,
and the non-competitive inhibition terms

φ 2 (x 2 ) = exp 1 -I in exp(-a(x 2 )h) K I -exp 1 -I in K I a(x 2 )h , ϕ(s 1 ) = K s K s + s 1 ,
where h is the pond depth and a(x 2 ) = a 0 + a 1 x 2 is the light extinction coefficient. It is straightforward to prove that this functions fulfill the assumptions H1-H3. The next table shows parameter values for the model defined by the system (3.6) based on [1]. 3.5 (from [1]) with a pond depth h = 0.1 m (left) and h = 0.5 m (right).

We consider the same benchmark cases as in [1]: an ammonium input concentration s in = 120 gN/m 3 , an incident light I in = 125 W/m 2 and two pond depths, either h = 0.1 m or h = 0.5 m.

For h = 0.1 m, we observe for low dilution rate (D < 0.17 d -1 ) that E a is globally stable. But for intermediate dilution rate (D ∈ (0.17; 1.1)), E a is unstable (given that λ s > λ n ) and the coexistence equilibrium E c is globally stable (see figure 3.5). For higher dilution rate, the nitrifiers are washed out and we obtain only algae (E a globally stable). For

h = 0.5 m, E c is globally stable if D < 1.1 d -1
, otherwise E a is globally stable (see figure 3.6). ) with initial condition (x 1 (0), x 2 (0), s 1 (0), s 2 (0)) = [START_REF] Berman | Nonnegative Matrices in the Mathematical Sciences[END_REF][START_REF] Gelhar | A critical review of data on field-scale dispersion in aquifers[END_REF][START_REF] Dutta | Fundamentals of biochemical engineering[END_REF][START_REF] Hsu | A survey of constructing Lyapunov functions for mathematical models in population biology[END_REF] and a pond depth h = 0.5 m, considering a dilution rate D = 0.9 (left) and D = 1.2 (right). In the left figure the equilibrium E c = (26.11, 95.34, 2.86, 109.5) is globally stable, however in the right figure the positive equilibrium E a = (0, 68.74, 114, 0) is globally stable.

Discussion

In this chapter we consider a model of two species, microalgae and nitrifying bacteria, in competition for nitrogen present as ammonium and nitrate produced by nitrification.

The system (3.6) defined in this work consider some important assumptions. First of all, we assume that the uptake functions of each specie are classical kinetics for substrate limitation (as Monod function for instance) but modified by assuming self shading and preference for ammo-nium over nitrate in the case of microalgae. With this we obtain density-dependent growth functions. Secondly, we assume cross-feeding, i.e., the nitrate produced in the reactor can be consumed by microalgae and we do not consider a initial input concentration of nitrate because is not reasonable in terms of the modelling.

We perform an analysis of existence and local stability of equilibria for the system (3.6), where conditions for existence and stability are given in terms of some parameters, in particular, dilution rate (D), ammonium input concentration (s in ) and break-even concentrations (see tables 3.2-3.3). To analyze the local stability of positive equilibrium was necessary to use a result of asymptotically autonomous systems to reduce dimension, obtaining the limiting system (3.18).

Concerning to global behaviour, the most interesting result was about ultimate bounded solutions near of coexistence equilibrium. This result guarantee that under some parameter conditions, the positive equilibrium can be globally asymptotically stable or at least the system solutions will remain close to this value. To prove this, it was necessary to use a two step approach. First, we realized that the system (3.18) can be interpreted as a perturbation of a two-species competition model for a single resource previously studied in [START_REF] Lobry | Effect on persistence of intra-specific competition in competition models[END_REF] where the authors provides a demonstration of global stability based on the Lyapunov direct method. Secondly, we can conclude the statement based on a non-vanishing perturbation result detailed in [START_REF] Khalil | Nonlinear Systems: Pearson New International Edition. Always learning[END_REF].

Finally, we shown an application of this model using parameter values from [1] where the authors conclude from simulations with a more complex model that nitrification does not occur for low dilution rate (D < 0.1 d -1 ) with h = 0.1 m. Our analysis is consistent with their results, although we obtain a different threshold.

Moreover, we can see that a deeper pond facilitates the presence of nitrifiers (E c becomes globally stable for all dilution rates lower than 1.1 d -1 ) given that it slows down the growth of microalgae (because of light attenuation). We can also observe that E n can never be locally stable within these operating conditions: microalgae can always invade a pure culture of nitrifiers.

Finally, this analysis also highlights the potential benefits of the presence of nitrifiers. For example, as it can be seen in Figure 3.4, the ammonium concentration with a dilution rate of 0.8 d -1 at the equilibrium with only algae (given by the green curves) is 75.3 gN/m 3 for h = 0.1 m and 111.3 gN/m 3 for h = 0.5 m (the concentration increases with pond depth given that algae growth becomes more light-limited). On the other hand, the ammonium concentration at the coexistence equilibrium (given by the blue curve) is 1.7 gN/m 3 , whatever is the pond depth. The presence of nitrifers allows a huge decrease of ammonium concentration at equilibrium. This is of particular interest in order to limit ammonia stripping (which can represent an undesirable loss of nitrogen if the objective is to produce microalgal biomass).

Our model and analysis do not include all the phenomena involved in the process (effect of O 2 , CO 2 , day-night cycle, ammonia stripping...). Nonetheless it gives us some hints on how to manage the presence of nitrifiers in an algae pond. This study also shows that a simple model can be useful to point out some phenomena from a mathematical analysis.

Appendices A Jacobian Matrix

The Jacobian matrix is the matrix of all first-order partial derivatives of a vector-valued function. Suppose f : R n → R m is a function (which takes as input real n-tuples and produces as output real m-tuples). Such a function is given by m real-valued component functions,

f 1 (x 1 , . . . , x n ), . . . , f m (x 1 , . . . , x n ).
The partial derivatives of all these functions with respect to the variables x 1 , . . . , x n can be organized in an m-by-n matrix, the Jacobian matrix J f of f , as follows:

J f (x 1 , . . . , x n ) =    ∂ ∂x 1 f 1 . . . ∂ ∂x n f n . . . . . . . . . ∂ ∂x 1 f m . . . ∂ ∂x n f m   
In the particular case when the function f : R n → R n define a dynamical system ẋ = f (x) with x ∈ R n , the Hartman-Grobman theorem [START_REF] Perko | Differential Equations and Dynamical Systems[END_REF] or Linearization theorem [6] assure that behavior of the system near a hyperbolic equilibrium point x * is the same to the linearized system ẋ = J f (x * )x.

For the system (3.6) the Jacobian matrix is

J(x 1 , x 2 , s 1 , s 2 ) =     J 11 0 J 13 0 0 J 22 J 23 J 24 -J 41 J 32 J 33 0 J 41 J 42 J 43 J 44     , (3.24)
where

J 11 = µ 1 (s 1 ) -D, J 13 = x 1 dµ 1 ds 1 (s 1 ), J 22 = ( ∂µ 2 ∂x 2 (s 1 , x 2 ) + ∂µ 3 ∂x 2 (s 1 , s 2 , x 2 ))x 2 + µ 2 (s 1 , x 2 ) + µ 3 (s 1 , s 2 , x 2 ) -D, J 23 = ( ∂µ 2 ∂s 1 (s 1 , x 2 ) + ∂µ 3 ∂s 1 (s 1 , s 2 , x 2 ))x 2 , J 24 = x 2 ∂µ 3 ∂s 2 (s 1 , s 2 , x 2 ), J 32 = -k 2 (x 2 ∂µ 2 ∂x 2 (s 1 , x 2 ) + µ 2 (s 1 , x 2 )), J 33 = -D -k 1 x 1 dµ 1 ds 1 (s 1 ) -k 2 x 2 ∂µ 2 ∂s 1 (s 1 , x 2 ), J 41 = k 1 µ 1 (s 1 ), J 42 = -k 2 (x 2 ∂µ 3 ∂x 2 (s 1 , s 2 , x 2 ) + µ 3 (s 1 , s 2 , x 2 )), J 43 = k 1 x 1 dµ 1 ds 1 (s 1 ) -k 2 x 2 ∂µ 3 ∂s 1 (s 1 , s 2 , x 2 ), J 44 = -D -k 2 x 2 ∂µ 3 ∂s 2 (s 1 , s 2 , x 2 ).

Chapter 4

Equivalence of finite dimensional input-output models of solute transport and diffusion

Introduction

In underground media (soil, unsaturated zone and aquifers), the transport of solutes is slow, remaining always slower than a few meters per hour in some large pore structures [START_REF] De Marsily | Quantitative hydrogeology: groundwater hydrology for engineers[END_REF][START_REF] Fetter | Contaminant Hydrogeology[END_REF]. Transport can also be orders of magnitude slower in narrow pores immediately neighboring the former larger pores giving rise to strong dispersive effects [START_REF] Gelhar | Stochastic Subsurface Hydrology[END_REF][START_REF] Gelhar | A critical review of data on field-scale dispersion in aquifers[END_REF]. Slow to very slow motions, high dispersion, and direct interactions between slow diffusion in small pores and fast advection in much larger pores are ubiquitous in soils and aquifers [START_REF] Coats | Dead-End Pore Volume and Dispersion in Porous Media[END_REF]. They are also the most characteristic features of underground transport as long as it remains conservative. Concerning reactivity, the characteristic feature of underground media is the strong water-rock interactions due both to the high surface to volume ratio and to the slow solute movements in the pores [START_REF] Steefel | Reactive transport modeling: An essential tool and a new research approach for the earth sciences[END_REF][START_REF] Van Genuchten | Mass transfer studies in sorbing porous media i. analytical solutions 1[END_REF]. The dominance of these characteristic features up to some meters to hundreds of meters have prompted the development of numerous simplified models starting from the double-porosity concept [START_REF] Warren | The behavior of naturally fractured reservoirs[END_REF]. In double-porosity models, solutes move quickly by advection in a first homogeneous porosity with a small volume representing focused fast-flow channels and slowly by diffusion in a second large homogeneous porosity. Exchanges between the two porosities is diffusion-like, i.e. directly proportional to the differences in concentrations. Such models have been widely extended to account not only for one diffusive-like zones but for many of them with different structures and connections to the advective zone [START_REF] Haggerty | Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity[END_REF][START_REF] Pruess | A Practical Method for Modeling Fluid and Heat Flow in Fractured Porous Media[END_REF]. Such extensions are thought to model both the widely varying transfer times and the rich water-rock interactions. The two most famous ones are the Multi-Rate Mass Transfer model (MRMT) [START_REF] Carrera | On matrix diffusion: formulations, solution methods and qualitative effects[END_REF][START_REF] Haggerty | Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity[END_REF] and Multiple INteracting Continua model (MINC) [START_REF] Narasimhan | An Approach for Analyzing Transport in Strongly Heterogeneous Systems[END_REF]. They are made of an infinity of diffusive zones deriving from analytic solutions of the diffusion equation in layered, cylindrical or spherical impervious inclusions (MRMT) or in series (MINC). Between the single and infinite diffusive porosities of the dual-porosity and these models, many intermediary models with finite numbers of diffusive porosities have been effectively used and calibrated on synthetic, field, or experimental data showing their relevance and usefulness [START_REF] Babey | Multi-Rate Mass Transfer (MRMT) models for general diffusive porosity structures[END_REF][START_REF] De Dreuzy | Influence of porosity structures on mixing-induced reactivity at chemical equilibrium in mobile/immobile multirate mass transfer MRMT and multiple interacting continua MINC models[END_REF][START_REF] Mckenna | Tracer tests in a fractured dolomite: 3. double-porosity, multiple-rate mass transfer processes in convergent flow tracer tests[END_REF][START_REF] Willmann | Coupling of mass transfer and reactive transport for nonlinear reactions in heterogeneous media[END_REF][START_REF] Zinn | Experimental visualization of solute transport and mass transfer processes in two-dimensional conductivity fields with connected regions of high conductivity[END_REF]. Theoretical grounds are however missing to identify classes of equivalent porosity structures, effective calibration capacity on accessible tracer test data, and influence of structure on conservative as well as chemically reactive transport. One can then naturally wonder which representation suits the best experimental data, and if the two particular MRMT and MINC models are not restrictive structures. This is exactly the problem we address in this work, from a theoretical approach based on linear algebra.

More precisely, we study the equivalence problem for a wide class of network structures and provide necessary and sufficient conditions, making the mathematical proofs explicit. We have a proof of the input-output equivalence of MINC and MRMT configurations in the simple case of three compartments (see Appendix A) but our aims of this analysis is to generalize this result for the case of n compartments. We stick to the framework of stationary flows (in the mobile zone) and assume water saturation in the immobile zones. Concretely, we consider a system of n compartments interconnected by diffusion, whose water volumes V i (i = 1 • • • n) are assumed to be constant over the time. One reservoir is subject to an advection of a solute. We shall called mobile zone this particular reservoir, and all the others n -1 reservoirs will be called immobile zones (see Fig. We aim at describing the time evolution of the concentrations S i (i = 1 • • • n) of the solute in the n tanks. The solute is injected in tank 1 with a water flow rate Q at a concentration S in , and withdrawn from the same tank 1 at the same water flow rate Q with a concentration S out = S 1 . Thus, the tank 1 plays the role of the mobile zone. We represent this system by a system of n iii. M ij ≤ 0 for any i = j iv. ∑ j M ij = 0 for any i

The diagonal terms of the matrix V represent the volumes of the n zones, and the off-diagonal terms M i,j of the matrix M are the (opposite) of the diffusive exchange rate parameters between zones i and j (equal to 0 if i is not directly connected to j). Properties i. and iv. are related to the connectivity of the graph between zones and the mass conservation (i.e. Kirchoff's law). One can proceed to the following reconstruction of matrices V and M from a given matrix A as follows.

Lemma 4.2 For a given matrix A that fulfills Assumptions 4.1, let π be a permutation of {1,

• • • , n} such that A π(i)π(i+1) = 0 (i = 1 • • • n) with π(1) = 1. Define the numbers V i as follows V π(i+1) = V π(i) A π(i)π(i+1) A π(i+1)π(i) i = 2 • • • n with V 1 = 1. Then M = -V(A + BB t )
where V is the diagonal matrix with V i as diagonal entries.

Matrices A that fulfill Assumption 4.1 are compartmental matrices, that have been extensively studied in the literature (see for instance [START_REF] Jacquez | Qualitative theory of compartmental systems[END_REF][START_REF] Walter | Compartmental Modeling with Networks[END_REF]). In the present work, we focus on properties for the specific structure of compartmental matrices that we consider.

Notations and preliminary results

For the sake of simplicity, we introduce the following notations

• for any vector X ∈ R n and matrix Z ∈ M n,n (R), we denote

X = [X i ] i=2•••n , Z = [Z ij ] i=2•••n j=2•••n .
• diag (X) denotes the diagonal matrix whose diagonal elements are the entries of the vector X

• we denote by

Vand (x 1 , • • • , v m ) the (square) Vandermonde matrix Vand (x 1 , • • • , v m ) =    1 x 1 • • • x m-1 1 . . . . . . . . . . . . 1 x m • • • x m-1 m    • we define the vectors in R n ✶ =    1 . . . 1    .
Lemma 4.3 Under Assumptions 4.1, the domain R n + is invariant by the dynamics for any non-negative control u. PROOF. Take a vector X that is on the boundary of R n + and set

I = {i ∈ 1 • • • n | X i = 0}
. At such a vector, one has Ẋi = ∑ j / ∈I A i,j X j + B i u , i ∈ I Note that the matrix A is Metzler (all its non-diagonal terms are non-negative) and B is a nonnegative vector. Consequently one has Ẋi ≥ 0 , i ∈ I which proves that any forward trajectory cannot leave the non-negative cone.

Remark Under Assumptions 4.1, the linear system (4.1) is positive in the sense that for any nonnegative initial state and non-negative control u(•), state and output are non-negative for any positive time (see [START_REF] Farina | Positive Linear Systems: Theory and Applications[END_REF]). PROOF. The matrix M is symmetric and consequently it is diagonalizable with real eigenvalues. Its diagonal terms are positive and off-diagonal negative or equal to zero. Furthermore one as

Mii = M i+1,i+1 = -∑ j =i+1 M i+1,j = -∑ j =i Mi,j -M i+1,1 ≥ -∑ j =i

Mi,j

The matrix M is thus (weakly) diagonally dominant. As each irreducible block of the matrix M has to be connected to the mobile zone (otherwise the matrix A will not be irreducible), we deduce that at least one line of each block has to be strictly diagonally dominant. Then each block is irreducibly diagonally dominant and thus invertible by Taussky Theorem (see [49, 6.2.27]). Finally the eigenvalues of the matrix M belong to the Gersgorin discs

G( M) = i λ ∈ R | |λ -Mi,i |≤ ∑ j =i | Mi,j |
and we deduce that each eigenvalues λi of M are positive. The matrix M is thus symmetric definite positive. Lemma 4.5 Under Assumptions 4.1, the matrix A is non singular. Furthermore, the dynamics admits the unique equilibrium ✶u, for any constant control u PROOF. Let X be a vector such that AX = 0. Then, one has BB t X = -V -1 MX or equivalently

MX = -V 1 X 1 B .
Definition 4.10 A matrix A that fulfills Assumptions 4.1 and such that the sub-matrix

à = [A ij ] i=2•••n j=2•••n
is diagonal is called a MRMT (multi-rate mass transfer) matrix.

MRMT matrices correspond to particular arrow structure of the matrix A:

A =               -Q V 1 -∑ i d 1i V 1 d 12 V 1 • • • • • • • • • d 1n V 1 d 12 V 2 -d 12 V 2 0 • • • • • • 0 . . . . . . . . . . . . d 1n V n 0 • • • • • • 0 -d 1n V n              
or star connections of depth one, where all the immobile zones are connected to the mobile one (see Fig. 4.2). In the context of general compartmental models, the flow digraph asociated to matrix A is an n-compartment mammillary system [2]. 00000000000000000000 00000000000000000000 00000000000000000000 00000000000000000000 00000000000000000000 00000000000000000000 11111111111111111111 11111111111111111111 11111111111111111111 11111111111111111111 11111111111111111111 11111111111111111111 MINC matrices correspond to particular structure:

S out V 1 V 2 V 3 V 4 V 5 d 12 S in Q Q V V 6
A =           Q V 1 -d 12 V 1 d 12 V 1 0 • • • • • • 0 d 12 V 2 -d 12 +d 23 V 2 d 23 V 2 0 • • • 0 . . . 0 • • • • • • 0 d n-1,n V n - d n-1,n V n          
As ✶u is an equilibrium of the system (4.1) for any constant control u, this Lemma allows then to claim the following result.

Lemma 4. [START_REF] Bertucco | Continuous microalgal cultivation in a laboratory-scale photobioreactor under seasonal day-night irradiation: experiments and simulation[END_REF] For any constant control u, ✶u is a globally exponentially stable of the dynamics (4.1).

Finally, we characterize the minimal MRMT representations as follows.

Lemma 4.14 For a minimal representation (A, B, C) where A is MRMT, the eigenvalues of the matrix à are distinct.

PROOF. The eigenvalues of à for the MRMT structure are

λ i = -d 1i /V i (i = 2 • • • n). If there exists i = j in {2, • • • , n} such that λ i = λ j = λ, one can consider the variable S ij = V i V i + V j S i + V j V i + V j S j
instead of S i , S j and write equivalently the dynamics in dimension n -1:

                             Ṡ1 = Q V 1 (S in -S 1 ) + ∑ k≥2,k =i,j d 1k V 1 (S k -S 1 ) + d 1ij V ij (S ij -S 1 ) . . . . . . Ṡk = d 1k V k (S 1 -S k ) k ∈ {2, • • • , n} \ {i, j} . . . . . . Ṡij = d 1ij V ij (S 1 -S ij ) with V ij = V i + V j and d 1ij = -(V i + V j )λ, which show that (A, B, C) is not minimal.
In the coming sections, we address the equivalence problem of any network structure that fulfills Assumption 4.1 with either a MRMT or a MINC structure. There are many known ways to diagonalize the sub-matrix à or tridiagonalize the whole matrix A to obtain matrices similar to A with an arrow or tridiagonal structure. The remarkable feature we prove is that there exist such transformations that preserve the signs of the entries of the matrices (i.e. Assumption 4.1 is also fulfilled in the new coordinates) so that the equivalent networks have a physical interpretation.

Equivalence with MRMT structure

We first give sufficient conditions to obtain the equivalence with MRMT. Proposition 4.15 Under Assumption 4.1, take an invertible matrix P such that P -1 ÃP = ∆, where ∆ is diagonal. If all the entries of the vector P -1 ✶ are non-null and the eigenvalues of à are distinct, the matrix R = 1 0 0 -P∆ -1 diag (P -1 A(2 : n, 1))

is invertible and such that R -1 AR is a MRMT matrix.

PROOF. Take a general matrix A that fulfills Assumption 4.1. From Lemma 4.6, Ã is diagonalizable with P such that P -1 ÃP = ∆ where ∆ is a diagonal matrix. Let G be the diagonal matrix

G = -∆ -1 diag (P -1 A(2 : n, 1))
and define R = PG. Note that one has A✶ = -B from Assumptions 4.1. The n -1 lines of this equality gives A(2 : n, 1) + Ã ✶ = 0 and one can write

P -1 A(2 : n, 1) + P -1 Ã ✶ = 0 ⇔ P -1 A(2 : n, 1) + P -1 ÃPP -1 ✶ = 0 ⇔ P -1 A(2 : n, 1) + ∆P -1 ✶ = 0
Thus having all the entries of the vector P -1 A(2 : n, 1) non-null is equivalent to have all the entries of the vector P -1 ✶ non null.

All the entries of the vector P -1 A(2 : n, 1) being non-null, R is invertible and one has

R-1 Ã R = G -1 P -1 ÃPG = G -1 ∆G = ∆.
One can then consider the matrix R ∈ M n,n defined as

R = 1 0 0 R with R -1 = 1 0 0 R-1 One has R -1 AR =   A 11 A(1, 2 : n) R R-1 A(2 : n, 1) ∆  
We show now that the matrix R -1 AR fulfills Assumptions 4.1.

One has straightforwardly

R -1 AR = -B t B -R -1 V -1 MR .
As the irreducibility of the matrix V -1 M is preserved by the change of coordinates given by X → R -1 X, Property i. is fulfilled.

The diagonal terms of -(R -1 AR + BB t ) are -A 11 -1 (which is positive) and the diagonal of -∆ which is also positive. Property ii. is thus satisfied.

We have now to prove that column R-1 A(2 : n, 1) and row A(1, 2 : n) R are positive to show Property iii. From the definition of the matrix G, one has

∆ = -G -1 diag (P -1 A(2 : n, 1)) = -diag ( R-1 A(2 : n, 1))
and thus one has R-1 A(2 : n, 1) = -∆ ✶

As the diagonal terms of ∆ are negative, we deduce that the vector R-1 A(2 : n, 1) is positive. As the matrix VA is symmetric, one can write V 11 A(1, 2 : n) = A(2 : n, 1) ′ Ṽ and then

V 11 A(1, 2 : n) R ′ = R′ ṼA(2 : n, 1) = -R′ Ṽ R∆ ✶
Note that the matrix R′ Ṽ R can be written T ′ T with T = Ṽ1/2 R, and that the matrix T diagonalizes the matrix S = Ṽ1/2 Ã Ṽ-1/2 :

T -1 ST = R-1 Ã R = ∆
The matrix S being symmetric, it is also diagonalizable with a unitary matrix U such that U ′ SU = ∆. As the eigenvalues of à are distinct, their eigenspaces are one-dimensional and consequently the columns of any matrix that diagonalizes S into ∆ have to be proportional to the corresponding eigenvectors. So the matrix T is of the form UD where D is a non-singular diagonal matrix. This implies that the matrix R′ Ṽ R is equal to D 2 , which is a positive diagonal matrix. As -∆ ✶ is a positive vector, we deduce that the entries of A(1, 2 : n) R are positive.

Note that ✶ is necessarily an eigenvector of R-1 (or R) for the eigenvalue 1: as one has A✶ = -B, one has also A(2 : n, 1) = -Ã ✶ and then

R-1 ✶ = -R-1 Ã-1 A(2 : n, 1) = -∆ -1 R-1 A(2 : n, 1) = -∆ -1 diag ( R-1 A(2 : n, 1)) ✶ = ✶ Finally, one has (R -1 AR + BB t )✶ = R -1 A✶ + B = -R -1 B + B = 0
which proves that Property iv. is verified.

We come back to the condition required by Proposition 4.15 and show that it is necessarily fulfilled for minimal representations (we recall from Lemma 4.9 that controllability implies a minimal representation in our framework). Proposition 4.16 Under Assumptions 4.1, the entries of the vector P -1 ✶ are non null for any P such that P -1 ÃP = ∆ with ∆ diagonal, when the pair (A, B) is controllable. Furthermore, the eigenvalues of à are distinct.

PROOF. From Lemma 4.6, Ã is diagonalizable with P such that P -1 ÃP = ∆ where ∆ is a diagonal matrix. Let X = P -1 ✶. One has

P -1 Ãk P = ∆ k =⇒ P -1 Ãk ✶ = ∆ k X , k = 1, • • • This implies P -1 ✶, Ã ✶, • • • , Ãn-1 ✶ = diag (X)Vand(λ 1 , • • • , λ n ) or equivalently P -1 C Ã, ✶ = diag (X)Vand(λ 1 , • • • , λ n-1 )
We deduce that when C Ã, ✶ is full rank, diag (X) and Vand(λ 1 , • • • , λ n-1 ) are non-singular, that is all the entries of X are non-zero and the eigenvalues λ 1 , • • • , λ n-1 are distinct. We show now that the controllability of the pair (A, B) implies that the pair ( Ã, ✶) is also controllable.

From the property A✶ = -B, one can write

A = A 11 L -Ã ✶ Ã
where L is a row vector of length n -1. Then one has

A✶ = -1 0 , A 2 ✶ = -A 11 Ã ✶ , A 3 ✶ = -A 2 11 + L Ã ✶ A 11 Ã ✶ + Ã2 ✶
that are of the form

A k ✶ = α k P k with P k = Ãk-1 ✶ + ∑ j≤k-2 β kj Ãj ✶ , for k = 2, 3
By recursion, one obtains

P k+1 = -α k à ✶ + Ãk ✶ + ∑ j≤k-2 β kj Ãj+1 ✶ = Ãk ✶ + ∑ j≤k-1 β k+1,j Ãj ✶ , for k = 2, • • • Then, one can write -C A,B = C A,A✶ = α 1 α 2 • • • α n 0 P 2 • • • P n from which one deduces rk (C A,B ) = n ⇒ rk (P 2 , • • • , P n ) = n -1 ⇒ rk ( Ã✶, • • • Ãn-1 ✶) = n -1
One can also write [ Ã✶, • • • Ãn-1 ✶] = ÃC Ã, ✶ and as à is invertible (Lemma 4.6), we finally obtain that C Ã, ✶ is full rank under controllability assumption.

Finally Propositions 4.15 and 4.16 lead to the following result.

Definition 4.19 (Lanczos algorithm) Let S be a symmetric matrix of size m and q 1 be a vector of norm equal to one. One defines the sequence π k = (β k , q k , r k ) as follows

• β 0 = 0, q 0 = 0, r 0 = q 1 ,

• if β k = 0, define q k+1 = r k /β k , α k+1 = q ′ k+1 Sq k+1 , r k+1 = (Sα k+1 I)q k+1β k q k and β k+1 = ||r k+1 ||.

One can straightforwardly check that the vectors q k provided by this algorithm are orthogonal and of norm equal to one. The algorithm stops for k < m ("breakdown") or k = m. A non-breakdown condition for this algorithm is given in [START_REF] Golub | Matrix Computations[END_REF]Th 10.1.1]: Proposition 4. [START_REF] Chisti | Biodisel from microalgae[END_REF] When rk (C S,q 1 ) = m, the sequence π k is defined up to k = m, and the matrix Q =

[q 1 • • • q m ] verifies Q ′ AQ =          α 1 β 1 0 β 1 . . . . . . . . . . . . . . . . . . . . . β m-1 0 β m-1 α m         
where the numbers

β i (i = 1 • • • m) are positive.
Lemma 4. [START_REF] Clarke | Nonsmooth Analysis and Control Theory[END_REF] The Lanczos algorithm applied to the matrix ∆ with q 1 = A(2 : n, 1)/||A(2 : n, 1)|| provides an orthogonal unitary matrix Q such that Q ′ ∆Q is symmetric tridiagonal with positive terms on the sub-(or super-) diagonal.

PROOF. As the matrix ∆ is diagonal, one has

C ∆,q 1 = Vand (λ 1 , • • • , λ n-1 )diag (q 1 )
where λ i (i = 1 • • • n -1) are the diagonal elements of ∆. Furthermore, as Assumptions 4.1 imply the equality A✶ = -B, one has

q 1 = - 1 ∑ n-1 i=1 λ i    λ 1 . . . λ n-1   
As λ i are all distinct and non null, q 1 is a non null vector and the matrices Vand (λ 1 , • • • , λ n-1 ), diag (q 1 ) are full rank. Therefore C ∆,q 1 is full rank and Proposition 4.20 can be used.

Let us recall the well known Cholesky decomposition of symmetric matrix.

Theorem 4.22

Let S be a symmetric definite positive matrix. Then, there exists an unique upper triangular matrix U with positive diagonal entries such that S = U ′ U.

We are ready now to explicit a tridiagonalization of the matrix A with positive entries on the sub-and super-diagonals. Proposition 4. [START_REF] Coats | Dead-End Pore Volume and Dispersion in Porous Media[END_REF] Let A be a MRMT matrix such that (A, B) is controllable. Let Q be the orthogonal matrix given by the Lanczos algorithm applied to ∆ with q 1 = A(2 : n, 1)/||A(2 : n, 1)||. Let U be the upper triangular matrix with positive diagonal entries given by the Cholesky decomposition of the symmetric matrix Q ′ ṼQ. Then the matrix

T = 1 0 0 QU -1
is such that T -1 AT is symmetric tridiagonal with positive entries on the sub-(or super-)diagonal.

PROOF. Lemma 4.21 provides the existence of the matrix Q such that Q ′ ∆Q is tridiagonal with positive terms on the sub-and super-diagonal. For convenience, we define the matrices

P = 1 0 0 Q and W = 1 0 0 U .
Clearly, P is orthogonal, W is upper triangular with positive diagonal, and one has T = PW -1 . Consider the matrix

P ′ AP = A 11 A(1, 2 : n)Q Q ′ A(2 : n, 1) Q ′ ∆Q .
For the particular choice of the first column of Q, one has

Q ′ A(2 : n, 1) = 1 ||A(2 : n, 1)||      1 0 . . . 0     
and Q ′ ∆Q is triangular with positive sub-diagonal. Therefore, P ′ AP is an upper Hessenberg matrix with positive entries on its sub-diagonal. Consider then

P ′ C A,B = P ′ B AB A 2 B • • • = P ′ B (P ′ AP)P ′ B (P ′ A 2 P)P ′ B • • •
Note that one has P ′ B = B and obtains recursively

P ′ B =         h 1 0 . . . . . . 0         , (P ′ AP)B =           ⋆ h 2 0 . . . . . . 0           , (P ′ A 2 P)B =          ⋆ ⋆ h 3 0 . . . 0          , • • •
where the number h i are positive. Therefore the matrix P ′ C A,B is upper triangular with positive diagonal entries, as the matrix W. Then T -1 C A,B = WP ′ C A,B is also upper triangular with positive entries on its diagonal. Proposition 6.22 (See part C in Annexes) implies that T -1 AT is tridiagonal with positive entries on its sub-diagonal. Let us show that T -1 AT is also symmetric. One has

T -1 AT = A 11 A(1, 2 : n)QU -1 UQ ′ A(2 : n, 1) UQ ′ ∆QU -1
As the matrix VA is symmetric by Assumption 4.1, one can write

A(1, 2 : n)QU -1 ′ = 1 V 1 (U -1 ) ′ Q ′ ṼA(2 : n, 1) = 1 V 1 (U -1 ) ′ U ′ UQ ′ A(2 : n, 1) = 1 V 1 UQ ′ A(2 : n, 1)
and as we have chosen V 1 = 1 we obtain A(1, 2 : n)QU -1 ′ = UQ ′ A(2 : n, 1). Consider now the sub-matrix UQ ′ ∆QU -1 . Note first that the decomposition Q ′ ṼQ = U ′ U implies the equalities

U ′ = Q ′ ṼQU -1 and (U -1 ) ′ = UQ ′ Ṽ-1 Q. Then on can write UQ ′ ∆QU -1 ′ = (U -1 ) ′ Q ′ ∆QU ′ = (UQ ′ Ṽ-1 Q)Q ′ ∆Q(Q ′ ṼQU -1 ) = UQ ′ Ṽ-1 ∆ ṼQU -1 = UQ ′ ∆QU -1
The matrix T provided by Proposition 4.23 possesses the following property.

Proposition 4. [START_REF] Coulson | Chemical and Biochemical Reactors and Process Control[END_REF] The vector X = T -1 ✶, where the matrix T is provided by Proposition 4.23, is positive.

PROOF. The matrices A + BB ′ and T -1 AT + BB ′ have non-negative entries outside their main diagonals. So there exists a number γ > 0 such that I + 1 γ (A + BB ′ ) and I + 1 γ (T -1 AT + BB ′ ) are non-negative matrices.

By Assumption 4.1, one has A✶ = -B, which implies the property

I + 1 γ (A + BB ′ ) ✶ = ✶ . Thus I + 1 γ (A + BB ′
) is a stochastic matrix (real square matrix, with each row summing to 1), and we know that its maximal eigenvalue is 1 (see for instance [START_REF] Berman | Nonnegative Matrices in the Mathematical Sciences[END_REF]Th 5.3]). As I + 1 γ (A + BB ′ ) and I + 1 γ (T -1 AT + BB ′ ) are similar:

T -1 I + 1 γ (T -1 AT + BB ′ ) T = I + 1 γ (A + BB ′ ), the maximal eigenvalue of I + 1 γ (T -1 AT + BB ′ ) is also 1. Furthermore, as A + BB ′ is irreducible by Assumption 4.1, I + 1 γ (T -1 AT + BB ′ ) is also irreducible. The property A✶ = -B implies I + 1 γ (T -1 AT + BB ′ ) X = X + 1 γ (T -1 A✶ + B) = X + 1 γ (-T -1 B + B) = X.
So X is an eigenvector of I + 1 γ (T -1 AT + BB ′ ) for its maximal eigenvalue 1. Finally note that X = T -1 ✶ implies that the first entry of X is equal to 1. Then by Perron-Frobenius Theorem (for nonnegative irreducible matrices, see for instance [START_REF] Berman | Nonnegative Matrices in the Mathematical Sciences[END_REF]Th 1.4])), we conclude that X is a positive vector.

We give now our main result concerning the MINC equivalence. PROOF. Let X = T -1 ✶ and Ā = R -1 AR. Define V = diag (X) 2 and M = -V( Ā + BB ′ ). As A + BB ′ is irreducible by Assumption 4.1, the similar matrix Ā + BB ′ is also irreducible, as well as M because V is a diagonal invertible matrix.

By Proposition 4.23, T -1 AT is a symmetric tridiagonal matrix with positive terms on the sub-or super-diagonal. By Proposition 4.24, X is a positive vector, and thus Ā = diag (X) -1 (T -1 AT)diag (X) is also a tridiagonal matrix with the same signs outside the diagonal. Thus, M is a tridiagonal matrix with negative terms on sub-or super-diagonal. Moreover, one has

M = -diag (X) 2 (diag (X) -1 T -1 ATdiag (X) + BB ′ ) = -diag (X)T -1 ATdiag (X) -X 2 1 .BB ′ .
where X 1 = 1. The matrix M is thus symmetric. One has

M✶ = -diag (X)T -1 ATdiag (X)✶ -BB ′ ✶ = -diag (X)T -1 ATX -B = -diag (X)T -1 AT(T -1 ✶) -B = -diag (X)T -1 A✶ -B = diag (X)T -1 B -B = diag (X)B -B = 0
The matrix M thus fulfills Assumption 4.1 and is tridiagonal: Ā is then a MINC matrix. Finally, one has B = R -1 B = B and C = CR = C.

Finally, Theorem 4.17 and Proposition 4.25 lead to the following result.

Examples: Reduction and the importance of minimal representation

Consider a network with one mobile zone and four immobile zones of identical volumes V i = 1 (i = 1 The structure of this network is neither MRMT nor MINC, and its corresponding matrix A is

u y V =1 2 1 3 4 V =1 V =1 V =1 5 V =1
A =       -3 1 1 0 0 1 -1 0 0 0 1 0 -3 1 1 0 0 1 -2 1 0 0 1 1 -2       .
One can easily compute the controllability matrix

C A,B =       1 -4 21 -129 906 0 1 -5 26 -155 0 2 -20 182 -1614 0 0 2 -18 136 0 0 6 -82 856      
and check that it is full rank (computing for instance det(C A,B ) = -896). Then the constructions of Sections 4.4 and 4.5 give the following equivalent MRMT and MINC matrices: We have checked numerically that each matrix A, A MRMT and A MI NC give the same co-prime transfer function T(z) = 14 + 47z + 45z 2 + 13z 3 + z 4 14 + 117z + 187z 2 + 92z 3 + 17z 4 + z 5 .

A MRMT =       - 4
Differently to the original network, the magnitude of the values of volumes and diffusive exchange rates are significantly different among compartments, opening the door of possible model reduction dropping some compartments. i. For the equivalent MRMT structure, one obtains ii. For the equivalent MINC structure, one obtains At the first look, this structure does not exhibit any special property or symmetry that could make believe that it is non minimal. By construction one has à ✶ = -✶ but the particular matrix A that we consider satisfies A(2 : 4, 1) = ✶. Consequently the vector A(2 : 4, 1) is an eigenvector of the matrix à for the eigenvalue -1. If the multiplicity of -1 was more than 1, then λ = -9.5 should be an eigenvalue of Ã, as the trace of à is -11.5. But an eigenvector X of à fulfills

V 1 = 1 , V 2 = 0.0398514 , V 3 = 0.0511169 , V 4 = 1 , V 5 =
V 1 = 1 , V 2 = 1.8 , V 3 = 1.
       Ṡ1 = -7S 1 + S 2 + 2S 3 + 3S 4 + u Ṡ2 = S 1 -7S 2 + 3S 3 + 3S 4 Ṡ3 = S 1 + 3 2 S 2 -5 2 S 3 Ṡ4 = S 1 + S 2 -2S 4 with the matrix A =     -7 1 2 3 1 -7 3 3 1 3 2 -5 2 0 1 1 0 -2     00000000
-7X 1 + 3X 2 + 3X 3 = λX 1 1.5X 1 -2.5X 2 = λX 2 X 1 -2X 3 = λX 3 one should have (λ + 7)X 1 -3X 2 -3X 3 = 0 with X 2 = 1.5 λ + 2.5 X 1 , X 3 = 1 λ + 2 X 1 (and X 1 = 0)
which is not possible for λ = -9.5. Then any matrix P that diagonalizes à should have one column proportional to the eigenvector õne, which amounts to have the vector P -1 ✶ with exactly one non-zero entry. Thus it is not possible to transform the system in a equivalent MRMT structure of the same dimension.

One can check that the pair (A, B) is indeed non controllable, even though the matrix à has distinct eigenvalues, as one has

AB =     -7 1 1 1     , A 2 B =     55 -8 -8 -8     = -B -8AB
from which one deduce rk (C A,B ) = 2. Indeed, the system admits a minimal representation of dimension 2 that can be found by gathering the immobile zones in one of volume V = V 2 + V 3 + V 4 = 6 and solute concentration

S = V 2 S 2 + V 3 S 3 + V 4 S 4 V = S 2 + 2S 3 + 3S 4 6
One can check that variables (S 1 , S) are solutions of the dynamics Ṡ1 = -7S

Discussion

We have shown that any general network structure is equivalent to a star structure (MRMT) or a series structure (MINC), that are commonly considered in geosciences to represent soil porosity in mass transfers. In this way, we reconcile these two different approaches, showing that they are indeed equivalent. Practically, this means that when the structure is unknown, or partially known, one can use equivalently the most convenient structure to identify the parameters or use some a priory knowledge.

Controllability property of a given mass transfer structure plays a crucial role. Although there is no particular control issue in the input-output representations of mass transfers, controllability is a necessary condition to obtain equivalence with the multi-rate mass transfer structures of depth one, introduced by Haggery and Gorelick in 1995 [START_REF] Haggerty | Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity[END_REF], or the multiple interacting continua structure. This condition is related to the minimal representation of linear systems, that is not necessarily fulfilled for such structures even for non-singular network matrices with distinct eigenvalues.

Although the objective of the present work is to show the exact equivalence of systems, we have shown on examples that MRMT and MINC representations could allow a simple and efficient way to obtain reduced models with a good approximation. Further investigations about such reduction techniques will the matter of a coming work.

From a geosciences view point, this analysis shows the existence of both identifiable and nonidentifiable porosity structures from input-output data in geological media where dispersion is controlled by diffusive exchanges between restricted but quick advective zones and slow but extensive diffusive zones. Input-output signals are typical of conservative tracer tests where nonreactive tracers are injected in an upstream well and analyzed in a downstream well [START_REF] Fetter | Contaminant Hydrogeology[END_REF]. Identifiable structures could thus be calibrated on tracer tests [2]. The porosity structure identified is however not unique as demonstrated on the example in Section 4.6, meaning that a porosity structure cannot be fully characterized by a conservative tracer test. While conservative transport and residence time distributions are already quite constraining for water-rock interactions [START_REF] De Dreuzy | Influence of porosity structures on mixing-induced reactivity at chemical equilibrium in mobile/immobile multirate mass transfer MRMT and multiple interacting continua MINC models[END_REF][START_REF] Donado | Multicomponent reactive transport in multicontinuum media[END_REF], highly non-linear reactivity depend on the structure of porosity beyond its input-output signature. In such cases, reactivity will not only depend on the input/output concentrations but also on the concentrations within the diffusion porosities, i.e. from the full state of the system. This should be seen as an advantage rather than as a drawback as some characteristic of the porosity structures might be revealed by appropriate combination of conservative and reactive tracers.

Appendices A Explicit equivalences for two models of three compartments

As a motivation of the general problem, we introduce the most simple case of configurations MINC and MRMT with three compartments (two inmobile and one mobile) and we will prove that this configurations are input-output equivalent using basic arguments of realization theory.

Let us begin with the MINC configuration that is presented in figure 4.8 One has the following dynamics for the concentrations S i with i = 1, 2, 3 in the three zones:

       . S 1 = Q V 1 (S in -S 1 ) + d 12 V 1 (S 2 -S 1 ) . S 2 = d 12 V 2 (S 1 -S 2 ) + d 23 V 2 (S 3 -S 2 ) . S 3 = d 23 V 3 (S 2 -S 3 ) (4.3) S S 1 in Q Q d 12 V 1 d V 2 V 3 S 3 S 2

23

S out We posit d = Q V 1 and

a 12 = d 12 V 1 , a 21 = d 12 V 2 , a 23 = d 23 V 2 , a 32 = d 23 V 3
Then the dynamic can be rewritten in the following way The transfer function of the input-output map S in → S out is given by

H C (s) = C(sI -A C ) -1 B.
In terms of the system parameters One has the following dynamics for the concentrations S i with i = 1, 2, 3 in the three zones:

H C (s) =
       . S 1 = Q V 1 (S in -S 1 ) + d 12 V 1 (S 2 -S 1 ) + d 13 V 1 (S 3 -S 1 ) . S 2 = d 12 V 2 (S 1 -S 2 ) . S 3 = d 23 V 3 (S 1 -S 3 ) (4.5) We posit d = Q V 1 and b 12 = d 12 V 1 , b 21 = d 12 V 2 , b 13 = d 13 V 1 , b 31 = d 13 V 3 .
Then the dynamic can be rewritten in the following way The transfer function of the input-output map S in → S out is given by

H C (s) = C(sI -A T ) -1 B.
In terms of the system parameters From the analysis above, we conclude that MINC and MRMT configurations with two inmobile compartments are input-output equivalent. Is it possible to generalize this result for the general case of n compartments? Naturally, it is not possible to replicate the previous process in the general case. For this reason, It was necessary to use another approach which was detailed in the previous sections.

H C (s) =

B Invariance of total volumes by MINC and MRMT transformations

Consider the minimal representation (A, B, C) with B and C as in system (4.1) assuming that A fullfill assumption 4.1. In the Theorem 4.17, using the transformation R, we obtain a algebraically Example For n = 4 in the system (4.1) consider the following matrices

A =     -2 1 0 0 3 -5 2 0 0 1 -2 1 0 0 2 -2     , B =     1 0 0 0     and C = 1 0 0 0 ,
where A is triadiagonal, fulfills Assumption H1 with

V MI NC =     1 0 0 0 0 1 3 0 0 0 0 2 3 0 0 0 0 1 3     and M =     1 -1 0 0 -1 5 3 -2 3 0 0 -2 3 4 3 -2 3 0 0 -2 3 2 3     ,
the system (A, B, C) are in MINC form. The controllability matrix of this system is

C 4 (A, B) =     1 0 0 0 -2 3 0 0 7 -21 3 0 -35 132 -27 6     ,
and det(C 4 (A, B)) = 54 = 0, so the system is controllable and the Theorem 4.17 ensures that there is a unique invertible transformation R such that R -1 AR have the MRMT form. By using a numerical algorithm and following the construction given by the Theorem, is possible find explicitly the transformation, in this case (the values of the coefficients are approximated) Example For n = 4 in the system (4.1) consider the following matrices

A =     -5 2 1 1 1 -1 0 0 3 0 -3 0 2 0 0 -2     , B =     1 0 0 0     and C = 1 0 0 0 ,
where A has a MRMT form and fulfills Assumption H1 with

V MRMT =     1 0 0 0 0 2 0 0 0 0 1 3 0 0 0 0 1 2     and M =     4 -2 -1 -1 -2 2 0 0 -1 0 1 0 -1 0 0 1     .
The controllability matrix of this system is

C 4 (A, B) =     1 -5 32 -210 0 1 -6 38 0 3 -24 168 0 2 -14 92     ,
and det(C 4 (A, B)) = 12 = 0, so the system is controllable and proposition 4.26 ensures that there is a unique invertible transformation R such that R -1 AR have the MINC form.

By using a numerical algorithm and following the construction given by the Theorem, is possible find explicitly the transformation, in this case (the values of the coefficients are approximated) An important observation from the examples above is that tr(V MI NC ) = tr(V MRMT ). This fact can have different interpretations depending on the problem under study. Suppose the MINC system is used to model the diffusion of the substance between several reactors with different volumes, then the MRMT model represents an equivalent model for the diffusion between the same amount of reactors with different volume distributions. However, both models maintain the same the total volume, from Proposition 4.27.

C A direct method to obtain an MINC structure

We consider now the system (A, B, C) defined in (4.1) under the same assumptions and notations, where the matrix A that fulfils Assumption 4.1.

The Multiple INteracting Continua models (MINC) matrix diffusion as diffusive-like exchanges within a succession of "continua" identified to the elementary cells issued by a finite-difference discretization of the diffusion process in the matrix. MINC has a chained-type connectivity where immobile zones are mutually connected on a line. The form of the matrix A MI NC is triadiagonal (not necessarily symmetrical).

The following result established a relation between the different models, where is possible find, under some assumptions, an invertible transformation such that A MI NC and A are similar and the tridiagonal matrix A MI NC fulfills Assumptions 4.1.

Theorem 4.28 Suppose that matrix A fulfills the assumption 4.1 and the realization (A, B, C) is controllable, then there exist a unique invertible matrix R such that A MI NC = R -1 AR fulfills assumptions 4.1 and is tridiagonal. PROOF. To construct the similarity transformation R we use the unsymmetric Lanczos algorithm. Let A ∈ M n such that A = V -1 (-BB t -M), A ii < 0, A ij ≥ 0 (i = j) and A✶ n = -B (the entries of every row of V -1 M sum to zero). Following the methodology of [START_REF] Golub | Direct block tridiagonalization of singleinput single-output systems[END_REF] the algorithm allows find a invertible square matrix P such that P -1 AP = H where H is a tridiagonal matrix, P -1 B = B and CP = C.

We consider the matrix C n (A, B) and through elementary row operations construct a matrix S 1 such that S 1 C n (A, B) = U 1 where U 1 is a upper triangular matrix with positive diagonal elements (algorithm 2). Note that this process is always possible, where C n (A, B) is invertible and from [42, Lemma 2.1] we have that S 1 AS -1

1 has a upper Hessenberg form and S 1 B = B and the elements under the principal diagonal are positives. Also, the S 1 matrix has the form

S 1 = 1 0 0 S1 ∈ M n,n ,
where S1 ∈ M n-1,n-1 is a lower triangular matrix. For simplicity, we dentote Ā = S 1 AS -1 1 .

Now, we consider the matrix C n ( Āt , B) that is invertible and through elementary row operations construct a matrix S 2 such that S 2 C n ( Āt , B) = U 2 where U 2 is a upper triangular matrix. Same as 108 above, we have that S 2 Āt S -1 2 has a upper Hessenberg form. Note that S 2 has a form similar to S 1 , CS t 2 = C and if we denote H t = S 2 Āt S -1 2 then H is a lower Hessenberg matrix.

If we define P = S -1 1 S t 2 by the construction above, H = P -1 AP is a tridiagonal matrix (upper and lower Hessenberg matrix), is non singular (similar to A invertible) and irreducible by construction, P -1 B = B, CP = C and the elements no nulls outside the principal diagonal are positives. Let R = PD a diagonal scaling (see remark C) where D = diag(P -1 ✶ n ) is a diagonal matrix and each diagonal entry is the row-sum of the matrix P -1 , we note that D✶ n = P -1 ✶ n ,

RB = R -1 B = B. Also, note that P -1 ✶ n = P -1 A -1 B = H -1 P -1 B = H -1 B,
is this, P -1 ✶ n is solution of the system HX = B and Lemma 6.21 above assure that each entry of the vector P -1 ✶ n is non null. Consecuently, each element of the principal diagonal of D is different to zero, ie D is invertible and we have that

R -1 AR = (PD) -1 A(PD) = D -1 (P -1 AP)D = D -1 HD,
where D and D -1 are diagonal matrices, R -1 AR is a tridiagonal matrix and the elements H i,i and H i+1,i H i,i+1 , i = 1, ..., n are invariant under this transformation (see remark C). Now

R -1 AR✶ n = R -1 APD✶ n = R -1 APP -1 ✶ n = R -1 A✶ n = R -1 (-B) = -B. Also R -1 AR = R -1 V -1 (-BB t -M)R = -R -1 V -1 BB t R -R -1 V -1 MR = -BB t -R -1 V -1 MR,
and R -1 V -1 MR is a tridiagonal matrix and fulfills the hyphotesis of the Proposition 4.27, then is possible find a matrix diagonal positive V such that V(R -1 V -1 MR) = M where M is a symmetric tridiagonal matrix. Summarizing

R -1 AR = -V-1 (BB t + M)
where M is irreducible symmetric tridiagonal matrix, Mii > 0, Mij ≤ 0, for i = j, M✶ n = 0. Finally A MI NC = R -1 AR is the matrix that that fulfils the statement.

Chapter 5

Conclusions and future perspectives

This chapter briefly reviews the main contributions of this thesis. The principal contribution of this thesis is the proposal and analysis of three different practical problems arising from real issues in bioprocess engineering. Each one is especially tailored and based on reasonable assumptions for each situation. The most specific results derived from this work can be divided into different topics depending of the studied problem.

Considering the first problem, chapter 2 was devoted to optimizing the productivity of microalgae in photobioreactors. The formulated model, stability results and optimization problem exhibit several features that can be of interest:

• In this batch model, we assume that there are no substrate limitations. Under this assumption, we focused on the influence of light in the growth of microalgae.

• We show that there is steady dynamical behavior over time on the part of the microalgae in the photobioractor, although in dark/ light cycles, the trajectories go to a quasi-periodic solution. This is not proved, but at least it is shown that there is a stable set where the solutions belong in time.

• It is feasible to find necessary optimality conditions for the optimization problem using Clarke generalized derivatives. Moreover, the results make sense under the model interpretations. For instance, proposition 2.9 indicates that the batch process should always conclude at the end of a period of light. That is, the respiration effects in the night period make the microalgae biomass decrease.

• We performed simulations using parameter values obtained from a preliminary study for microalgae chlamydomonas reinhardtii , and the results obtained are consistent.

These features and the methodology employed, which is based on ordinary differential equations with a discontinuous right side, nonsmooth optimization and numerical analysis, appear to be well suited to the scope of a thesis project.

The second problem was studied in chapter 3 and addressed modelling and stability analysis of a microalgal pond with nitrification. The main novelties and results derived from the analysis are given below.

• In this chemostat model, an intra-specific competition phenomenon was considered based on density-dependent growth functions and cross-feeding, i.e., the nitrate produced by the nitrifiers can be consumed by the algae. The last assumption is not generally considered in competition models.

• The coexistence of species in the chemostat was shown, although their stability may vary depending on certain parameter values. A Lyapunov stability approach was necessary to obtain this result.

• It is feasible to reduce the system using the theory of asymptotically autonomous systems, and the limiting system obtained can be considered a perturbation of a system of two species competing for a substrate. Using a result of non-vanishing perturbated systems was possible to obtain a strong stability theorem for equilibrium coexistence.

• We performed simulations using functions and parameter values obtained from a preliminary study [1], and the results obtained are equivalents from a qualitative point of view, but our model is much simpler than that used by the authors, where the referred work is a purely numerical study and based on simulations.

We take into account the third problem, which was considered in chapter 4 and address the equivalence of finite dimensional input-output models of solute transport and diffusion; the main conclusions are as follows:

• We have shown that any general network structure is equivalent to a star structure (MRMT) or a series structure (MINC), which are commonly considered in geoscience problems to represent soil porosity in mass transfers.

• The controllability property of a given mass transfer structure plays a crucial role because it is a necessary condition to prove observability in the special matrices of linear compartmental systems that are involved in the modelling. As a consequence, the minimality of the transfer function associated with the system implies the uniqueness of equivalent transformation.

• Examples show that MRMT and MINC representations could allow a simple and efficient way to obtain reduced models with good approximation.

The contributions of this thesis create possibilities for future work. Further from the perspectives discussed here, which originate from our contributions, a number of related research directions are worth exploring. This chapter briefly discusses the most important ones. Some of these research proposals represent on going work.

For the first problem, some possible perspectives are the following:

• Decrease the number of assumptions. For instance, by avoiding the assumption about substrate limitations, the growth model will have a larger dimension considering the substrate concentration equation.

• Consider different assumptions about natural light incidence. What happens if we consider a Fourier series approach instead of a step function approach to model superficial light in the photobioreactor process?

• What would happen in the case of species competition in the photobioreactor in a batch or continuous process?

For the second problem, some possible perspectives are as follows:

• Propose an optimal control problem to maximize the productivity of the microalgae subject in the system under consideration, where the dilution rate is the control variable, in order to obtain an optimal dilution profile for process management.

• Consider microalgae growth limitations as result of excess ammonia, which is more realistic and has not been considered in this thesis.

For the third problem, some possible perspectives are as follows:

• Consider input and output flows with different diffusion values in each compartment. What happens in the MIMO case?.

• Analyze the MINC and MRMT gradostat with three bioreactors and biomass diffusion.

• Generalize the result of the last analysis to n vessels, considering biomass diffusion with different growth functions (linear, Monod, Haldane, etc.) as compared to the case without biomass diffusion. Is there a connection between systems with and without biomass dynamics?

To conclude, we would like to stress a fact further elicited by this work, namely, that mathematical models in bioprocess modelling provide theoretical results that offer an idea of the process behavior and this helps in management decision making. We believe that the application of similar approaches may prove productive for problems that arise in other bioprocess issues. Suppose we have a linear system with transfer function H(s) (which can be a matrix, in general). We have seen that the transfer function is related to the matrices in the state space model via (6.2). Recall that the transfer function describes how the input to the system affects the output (when the initial state of the system is zero). In some sense, this might seem to indicate that the exact representation of the internal states of the system might be irrelevant, as long as the inputoutput behavior is preserved. We will see that there are multiple state-space realizations for a given system that correspond to the same transfer function. Consider any particular state-space model of the form (6.1). Now, let us choose an arbitrary invertible n × n matrix T, and define a new state vector x(t) = Tx(t)

In other words, the states in the vector x(t) are linear combinations of the states in the vector x(t). Since T is a constant matrix, we have Thus the transfer function for the realization with state-vector x is the same as the transfer function for the realization with state-vector x. For this reason, the transformation x = Tx is called a similarity transformation and systems say (algebraically) equivalent. Since T can be any invertible matrix, and since there are an infinite number of invertible n × n matrices to choose from, we see that there are an infinite number of realizations for any given transfer function H(s). From input-output viewpoint, two equivalent dynamical equations give the same transfer function and this remain unchanged under an similartity transformation. In other words, algebraically equivalent realizations have identical controllability, observability and stability properties.

The system defined in (6.1) can be writing In general, all linear system can be represented by

P(∆)x = Q(∆)u (6.3)
where x ∈ R n , u ∈ R m , ∆ is the differentiation operator, P(∆) and Q(∆) are, respectively, n × n and n × m matrices whose entries are polynomials in operator ∆. We assume that det(∆) = 0. Let us assume that the maximal degree of the entries of P is d p and the maximal degree of the entries of Q is d q . Let x (k) (0) = 0 for k = 0, ..., d p and u (k) (0) = 0 for k = 0, ..., d q . Applying the Laplace transform on both sides of the input-output equation (6.3) we obtain P(s)X(s) = Q(s)U(s), where X(s) = L{x}(s) and U(s) = L{u}(s) are the Laplace transforms of x and u, respectively. As P(s) is invertible for almost all s, we get X(s) = P(s) -1 Q(s)U(s). As usually, H(s) = P(s) -1 Q(s) is called the transfer matrix of system.

In [3,[START_REF] Chen | Linear System Theory and Design[END_REF], the author give the next definition

C The unsymmetric Lanczos procedure and tridiagonalization

The Lanczos biorthogonalization algorithm is an extension to nonsymmetric matrices of the symmetric Lanczos algorithm. The nonsymmetric Lanczos algorithm is quite different in concept from Arnoldiâ Ȃ Źs method because it relies on biorthogonal sequences instead of orthogonal sequences. The algorithm proposed by Lanczos for nonsymmetric matrices builds a pair of biorthogonal bases for the two subspaces K m (A, v 1 ) = span{v 1 , Av 1 , • • • , A n-1 v 1 } and K m (A t , w 1 ) = span{w 1 , A t w 1 , • • • , (A t ) n-1 w 1 }.

The algorithm that achieves is the following Algorithm 1 The Lanczos Biorthogonalization procedure Choose two vectors v 1 , w 1 such that (v 1 , w 1 ) = 1, Set β 1 = δ 1 = 0, w 0 = v 0 = 0, For j = 1, 2, ..., m Do: α j = (Av j , w j ), vj+1 = Av jαv jβ j v j-1 , ŵj+1 = A t w jαw jδ j w j-1 ,

δ j+1 = |( vj+1 , ŵj+1 )| 1 2
, If δ j+1 = 0 Stop β j+1 = ( vj+1 , ŵj+1 )/δ j+1 , w j+1 = ŵj+1 /β j+1 , v j+1 = vj+1 /δ j+1 End Do Note that there are numerous ways to choose the scalars δ j+1 , β j+1 . These two parameters are scaling factors for the two vectors v j+1 and w j+1 and can be selected in any manner to ensure that (v j+1 , w j+1 ) = 1. As a result of the algorithm, it is only necessary to choose two scalars δ j+1 , β j+1 that satisfy the equality δ j+1 β j+1 = ( vj+1 , ŵj+1 ). (

The choice taken in the above algorithm scales the two vectors so that they are divided by two scalars which have the same modulus. Both vectors can also be scaled by their 2-norms. In that case, the inner product of v j+1 and w j+1 is no longer equal to 1 and the algorithm must be modified accordingly. Consider the case where the pair of scalars δ j+1 , β j+1 is any pair thart satisfies the relation (6.5). Denote by T m the tridiagonal matrix

T m =       α 1 β 2 δ 2 α 2 β 3 . . . δ m-1 α m-1 β m δ m α m       We have                    a 1 b 1 0 • • •
and the equation c i-1 x i-1 + b i x i+1 = 0. (6.9)

We note that (6.7)-(6.9) is a overdetermined system, but have a unique solution (from the hypothesis) and this implies that the solution of the system (6.8) necessarily is not the trivial, ie exist two linearly dependent rows and where the system (6.8) is homogeneous, then is the trivial solution plus an infinite set of other solutions. This is a contradiction where the solution of the system is unique. Therefore, each entrie x i , (i = 1, . . . , n) of the vector solution X is different to zero. This complete the proof.

We also recall a nice result about tridiagonalization of single-input single-output systems, from [START_REF] Golub | Direct block tridiagonalization of singleinput single-output systems[END_REF].

DTD -1 =          a 1 √ b 1 c 1 0 • • • 0 √ b 1 c 1 a 2 √ b
• • • 0 b n-1 c n-1 a n         
Since DTD -1 is symmetric and T and DTD -1 have the same eigenvalues, this implies that the eigenvalues of T are real. Suppose λ is an eigenvalue of T. Then λI -T is also of the same form that the matrix T. If the first row and last column of λI -T are deleted, then the resulting matrix is an n -1 × n -1 upper triangular matrix with no zero entries on its main diagonal, since b i c i > 0, for i = 1, 2, ..., n -1. Hence this submatrix has rank n -1. It follows that λI -T has rank at least n -1. However, λI -T has rank at most n -1 since λ is an eigenvalue of T. So by definition λ has geometric multiplicity one.

D Compartmental Systems

Kinetics is that branch of dynamics that pertains to the turnover of specific particles in a biological system. A compartment is an amount of material which acts kinetically in a homogeneous distinct way. The compartment to which a particle belongs characterizes both its physical-chemical properties and its environment. A compartment may not be an actual physical volume; however, the amount of some material in a physiological space is often treated as a compartment in certain clinical studies. The particles of each compartment are influenced by forces which cause the particles to transfer from one compartment to another. All particles in a particular compartment have the same probability of transition since within the compartment all particles are well-mixed and considered indistinguishable by the system. The transition from one compartment to another occurs by passing through some physical barrier or by undergoing some physical or chemical transformation. Associated with each compartment is its size. Sometimes the terms volume and size are used interchangeably.

A compartmental system consists of two or more compartments, interconnected in the sense that among certain compartments there is exchange of material. The compartmental system will be primarily modeled in a continuous deterministic manner by a collection of ordinary differential equations, each equation describing the time rate of change of amount of material in a particular compartment. These rates of change are dictated by the physical-chemical laws that govern the material exchange between interacting compartments, e.g., diffusion, temperature, chemical reactions, etc. It will be assumed that the set of equations model average properties of the system for very large numbers of particles. From the idealized physical model it is hypothesized as to the exact nature of the interconnections that occur between compartments. In the usual representation for a compartmental system, a box denotes a compartment, and an arrow indicates the transfer of material into or out of a compartment. There also may be inputs from the outside environment into one or more compartments (vertical arrows pointing into the tops of boxes), and there can be excretion of material from some of the compartments to the outside environment (vertical arrows pointing out of the bottoms of boxes). Should there be no exchange of material to the outside environment, the compartmental system is referred to as closed; otherwise it is said to be an open system. Realistically, many compartmental systems are open, for some material is lost by excretion, metabolism, etc.

Linear compartmental models

The kinetics of the tracer amount x i (t) will be assumed to follow the constant coefficient linear equation ẋi (t) = b i (t) + n ∑ j=1,j =i a ij x j -n ∑ j=0,j =i a ji x i , j = 1, ..., n.

(6.10)

The b i (t) is the input rate of tracer to compartment i from outside the system. The constant fractional transfer coefficients a ij (j = i) are nonnegative. The last sum in (6.10) can be rearranged There are some important types of connectivity in compartmental models that appear repeatedly in idealized physiological systems. The first is the catenary system in which the compartments connect in a series or chain as is shown in figure 6.2. Only 'adjacent compartments communicate. Of special interest are catenary systems wherein there is tracer input into either compartment 1 or n, and single exit only from one of those two compartments. In particular, a common physiological catenary 3-compartment model is with the compartments identified as (1) blood plasma; (2) interstitial fluid; (3) cells. Here, due primarily to the kidney, the first compartment has the single exit to the outside. The compartmental matrix A for a catenary system as shown in figure 6.2 is easily seen to be tridiagonal. Often useful in analyzing the kinetics of the distribution of a tracer injected into blood plasma and which enters the interstitial space is the mammillary system (see figure 6.3). Here one compartment acts as the "mother" or central compartment and all other compartments are "daughters". Connectivity only takes place between the mother compartment and each individual daughter. It is standard to number the compartments commencing with the central compartment designated as number one. Provided this numbering scheme is instituted, the compartmental matrix A has nonzero entries only in the first row, first column, and on the main diagonal. A complete analysis about solution of this systems and stability, strongly connected property and irreducibility, non negative eigenvalues of compartmental matrices is given in [2]. These results will be useful in Chapter 4 of this thesis.

2 . 1

 21 Parameters and functions used in the general model. . . . . . . . . . . . . . . . . . . 2.2 optimal values in different light environments for C. reinhardtii considering constant light. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Parameter estimations for different incidental light settings for C. reinhardtii. . . . . . 2.4 Optimal values in different incidental light settings for C. reinhardtii considering constant light. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Parameter values for C. reinhardtii obtained from [58]. . . . . . . . . . . . . . . . . . . 2.6 Linear regression f (x) = axL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Parameters estimated for C. reinhardtii in this model. . . . . . . . . . . . . . . . . . . .3.1 Variables used in the algal pond model . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Equilibrium existence conditions for break-even concentrations . . . . . . . . . . . . 3.3 Local stability conditions of non-coexistence equilibrium points for break-even concentrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Equilibrium existence and local stability conditions for break-even concentrations in system (3.18). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Parameter values for the system (3.6). . . . . . . . . . . . . . . . . . . . . . . . . . . . xv List of Figures 1.1 Effect of nutrient concentration on the specific growth rate of E. coli (source [87]). . . 1.2 Batch bioreactor in a laboratory experiment (source Laboratoire de Biotechnologie de l'Environnement, INRA, Narbonne, France ). . . . . . . . . . . . . . . . . . . . . . 1.3 Growth of phytoplankton in a continuous reactor (source Laboratoire d'Ecologie du plancton, Marin, CNRS, Villefranche sur mer, France). . . . . . . . . . . . . . . . . . . 1.4 Batch growth profile (source [17]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Light-response curve of photosynthesis (P-curve). The intercept on the vertical axis is the measure of O2 uptake due to dark respiration. I c , light compensation point; I s , light saturation intensity; I h , light intensity value at which photoinhibition occurs (source [84]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 An open raceway pond is a photobioreactor with natural light (source Laboratoire de Biotechnologie de l'Environnement, INRA, Narbonne, France ). . . . . . . . . . . 1.7 Multiple chemostats connected in a laboratory experiment (source Laboratoire de Biotechnologie de l'Environnement, INRA, Narbonne, France ). . . . . . . . . . . . . 1.

2 (

 2 right). In the left figure the equilibrium E c = (26.11, 95.34, 2.86, 109.5) is globally stable, however in the right figure the positive equilibrium E a = (0, 68.74, 114, 0) is globally stable. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Example of a network with one mobile zone . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Example of a MRMT network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Example of a MINC network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Example of a network with one mobile and four immobile zones . . . . . . . . . . . 4.5 Nyquist diagrams (black: original system, blue: reduced MRMT, green: reduced MINC) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Structure of the example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Simplified equivalent structure of the example . . . . . . . . . . . . . . . . . . . . . . 4.8 The MINC configuration with two inmobile compartments. . . . . . . . . . . . . . . 4.9 The MRMT configuration with two inmobile compartments. . . . . . . . . . . . . . . 6.1 The Compartmental System Connectivity Diagram (source [2]). . . . . . . . . . . . . 6.2 An n-Compartment Catenary System (source [2]). . . . . . . . . . . . . . . . . . . . . 6.3 An n-Compartment Mammilary System (source [2]). . . . . . . . . . . . . . . . . . . . xvii xviii Chapter 1

Figure 1 . 1 .

 11 Figure 1.1. Effect of nutrient concentration on the specific growth rate of E. coli (source [87]).

Figure 1 . 2 . 0 -

 120 Figure 1.2. Batch bioreactor in a laboratory experiment (source Laboratoire de Biotechnologie de l'Environnement, INRA, Narbonne, France ).

Figure 1 . 3 .

 13 Figure 1.3. Growth of phytoplankton in a continuous reactor (source Laboratoire d'Ecologie du plancton, Marin, CNRS, Villefranche sur mer, France).

) 3 . 4 .

 34 Mixed inhibition: both equations (1.16) and (1.17) apply. Pure noncompetitive inhibition: mixed inhibition with K IC = K IU . where K app A represents the apparent value of specificity coefficient K A [L mol -1 s -1 ], I represents the inhibitor concentration [mol L -1 ], K IC represents the competitive inhibition constant for I [mol L -1 ], K IU represents the uncompetitive inhibition constant for I [mol L -1 ] and K app 0 represents the apparent value of catalytic constant K 0 [s -1 ].

A

  . Adaptation (lag phase) B. Exponential growth phase (log phase) C. Stationary phase D. Logarithmic death phase

Figure 1 . 4 .

 14 Figure 1.4. Batch growth profile (source [17]).

Figure 1 . 5 .

 15 Figure 1.5. Light-response curve of photosynthesis (P-curve). The intercept on the vertical axis is the measure of O2 uptake due to dark respiration. I c , light compensation point; I s , light saturation intensity; I h , light intensity value at which photoinhibition occurs (source [84]).

Figure 1 . 6 .

 16 Figure 1.6. An open raceway pond is a photobioreactor with natural light (source Laboratoire de Biotechnologie de l'Environnement, INRA, Narbonne, France ).

Figure 1 . 7 .

 17 Figure 1.7. Multiple chemostats connected in a laboratory experiment (source Laboratoire de Biotechnologie de l'Environnement, INRA, Narbonne, France ).

Figure 1 . 8 .

 18 Figure 1.8. The standard n-vessel gradostat. The left vessel labeled R is a reservoir containing nutrient at concentration S 0 , C is an overflow vessel, and D denotes the dilution rate. All vessels have the same volume (source [88]).

Figure 1 . 9 .

 19 Figure1.9. Irreducible "dead-end" gradostat. Note that the inflow to each vessel balances the outflow (source[START_REF] Smith | The Theory of the Chemostat[END_REF]).

Figure 2 . 1 .

 21 Figure 2.1. Monod growth function.

Figure 2 . 2 .

 22 Figure 2.2. Phase plane of the differential equation (2.7) for microalgae C. reinhardtii with parameter values µ = 2.34 [d -1 ], ρ = 0.34 [d -1 ], c = 0.253 [g.L -1 ]. and different initial conditions.

Figure 2 . 3 .

 23 Figure 2.3. Solutions of the differential equation (2.9) for microalgae C. reinhardtii with parameter values µ = 2.34 [d -1 ], ρ = 0.34 [d -1 ], c = 0.253 [g.L -1 ]. considering regular intervals of 12 hours (0.5 [d]) and different initial conditions. A stable set exist in the interior of the interval [2.35, 2.8] over biomass concentration axis.

Figure 2 . 4 .

 24 Figure 2.4. Surface and level curves of the net rate of production (mean biomass volumetric productivity) for the optimization problem (2.6) for C. reinhardtii at the parameter values µ = 2.34 [d -1 ], ρ = 0.34 [d -1 ], c = 0.253 [g.L -1 ].

Figure 2 . 5 .

 25 Figure 2.5. Surface and level curves of the net rate of production for the optimization problem (2.6) for C. reinhardtii with the same parameter values in summer period.

Figure 2 . 6 .

 26 Figure 2.6. Surface and level curves of the net rate of production for the optimization problem (2.6) for C. reinhardtii with the same parameter values in regular time intervals.

Figure 2 . 7 .

 27 Figure 2.7. Surface and level curves of the net rate of production for the optimization problem (2.6) for C. reinhardtii with the same parameter values in winter period.

Figure 2 . 8 .

 28 Figure 2.8. Biomass trajectories associated to the optimal values of initial concentration x * 0 and terminal time T * in different light environments, with turnaround times t a = 1 (left) and t a = 2 (right).

  Figure 2.9. The figure shows the results of the table2.4. It may be noted that the higher the incidencial light, higher the final biomass concentration, while the terminal batch time decreases and the mean volumetric productivity increases, the initial concentration needed to achieve optimal productivity decreases slightly.

I 0

 0 t)+K I I 0 (t)e -axL +K I s s+K s , Replacing all previously considered in the model (2.1) we obtain the system  (t)+K I I 0 (t)e -axL +K I s s+K s xρx, (2.13)

  terms of the unit measurement used in this model, the last value is equivalent to ρ = 0.014 [h -1 ].

Figure 2 . 10 .

 210 Figure 2.10. Linear regression and zero intercept for the estimation of the light attenuation parameter a.

Figure 2 . 11 .

 211 Figure 2.11. Approximation of functions µ 1 (.) and µ 2 (.) by nonlinear least squares method.

Figure 3 . 1 .

 31 Figure 3.1.In this setting μn = 2, μs = 1, K n = 1, K s = 0.12, K x = 0.5, s in = 1, D = 0.5, k 1 = 1 and considering the growth functions in the application section, the graph show the existence of the three equilibrium points that was proven in lemma above.

Figure 3 . 2 .

 32 Figure 3.2. Graph of the function µ 3 (.) for k 2 = 1 and the parameter values in table 3.5.

Figure 3 . 4 .

 34 Figure 3.4. Growth functions with parameter values of table3.5 (from[1]) with a pond depth h = 0.1 m (left) and h = 0.5 m (right).

3 )Figure 3 . 5 .

 335 Figure 3.5. Solutions of the system (3.6) with initial condition (x 1 (0), x 2 (0), s 1 (0), s 2 (0)) = (10, 40, 30, 50) and a pond depth h = 0.1 m, considering a dilution rate D = 0.1 (left) and D = 0.4 (right). In the left figure the equilibrium E a = (0, 1496, 0.045, 0) is globally stable, however in the right figure the positive equilibrium E c = (11.2461, 1034, 0.388, 36.88) is globally stable.

3 )Figure 3 . 6 .

 336 Figure 3.6. Solutions of the system (3.6) with initial condition (x 1 (0), x 2 (0), s 1 (0), s 2 (0)) = (10, 40, 30, 50) and a pond depth h = 0.5 m, considering a dilution rate D = 0.9 (left) and D = 1.2 (right). In the left figure the equilibrium E c = (26.11, 95.34, 2.86, 109.5) is globally stable, however in the right figure the positive equilibrium E a = (0, 68.74, 114, 0) is globally stable.

Figure 4 . 1 .

 41 Figure 4.1. Example of a network with one mobile zone

Lemma 4 . 4

 44 Under Assumptions 4.1, the matrix M is symmetric definite positive.

Figure 4 . 2 .

 42 Figure 4.2. Example of a MRMT network

Proposition 4 .

 4 [START_REF] De Dreuzy | Influence of porosity structures on mixing-induced reactivity at chemical equilibrium in mobile/immobile multirate mass transfer MRMT and multiple interacting continua MINC models[END_REF] Let A be a MRMT matrix such that (A, B) is controllable and R = Tdiag (T -1 ✶), where T is provided by Proposition 4.23. Then (R -1 AR, B, C) is an equivalent representation where R -1 AR is a MINC matrix.

Theorem 4 . 26

 426 Any minimal representation (A, B, C) that fulfills Assumptions 4.1 is equivalent to a MINC structure. Remark From Theorem 4.26 and following the same argument that in Corollary 4.18 for system (4.1) the minimality of the realizations (A, B, C) and (A MI NC , B, C) implies that this realizations are input-output equivalents.

Figure 4 . 4 .

 44 Figure 4.4. Example of a network with one mobile and four immobile zones

  2.9090317 with d 12 = 0.3256267 , d 13 = 0.1692779 , d 14 = 1, d 15 = 1.5050954and notes that zones 2 and 3 are of relatively small volumes (compared to the total volume of the system which is equal to 5) and connected to the mobile zone with relatively small diffusive parameters. Then one may expect to have a good approximation with a reduced MRMT model dropping zones 2 and 3. Keeping the volumes V 1 , V 4 , V 5 with the parameters d 14 , d 15 , one obtains the 3 compartments MRMT matrix + 1.5173871z + z 2 0.5173871 + 4.8359736z + 5.0224825z 2 + z 3 .

  6666667 , V 4 = 0.3603604 , V 5 = 0.1729730 with d 12 = 3 , d 23 = 6 , d 34 = 0.8888889, d 45 = 0.1636231 .Here, one notes that the two last volumes are relatively small and connected with relatively small diffusion terms. Keeping the volumes V 1 , V 2 , V 3 with the parameters d 12 , d 23 , one obtains the

Figure 4 . 6 .

 46 Figure 4.6. Structure of the example

Figure 4 . 8 .

 48 Figure 4.8. The MINC configuration with two inmobile compartments.

  S + BS in S out = CS with S = (S 1 , S 2 , S 3 ) T and

(4. 4 )Figure 4 . 9 .

 449 Figure 4.9. The MRMT configuration with two inmobile compartments.

Remark

  S + BS in S out = CS with S = (S 1 , S 2 , S 3 ) T and A T =   -d -b 12 -b 13 b 12 Configurations (d, b 12 , b 13 , b 21 , b 23 ) and d, b12 , b13 , b21 , b23 with d = d, b21 = b 13 , b13 = b 12 , b21 = b 31 and b31 = b 21 are equivalents (one has just to interchange subindices 2 and 3).

ds 2 31 ( 4 . 6 )

 23146 +d(b 21 +b 31 )s+db 21 b 31 s 3 +(b 12 +b 13 +b 21 +b 31 +d)s 2 +(b 12 b 31 +b 21 b 13 +b 21 b 31 +d(b 21 +b 31 ))s+db 21 b One can check that the pairs (A T , B) and (A T , C) are respectively controllable and observable exactly when b 21 = 0, b 31 = 0 and b 21 = b 31 . Under these conditions, numerator and denominator of H T are coprimes and the realization is thus minimal.

  (A MRMT , B, C) is controllable because both minimal realizations are similar.

  (A MI NC , B, C) is controllable because both minimal realizations are similar as in the anterior example.

••

  The observability matrix associated to the pair (A, C) is given by O A,C = A system Ẋ = AX + Bu is said to be controllable when rk (C A,B ) = n, and observable for y = CX when rk (O A,C ) = n.

  ) = TAx(t) + BTu(t) = TAT -1 x(t) + TBu(t)y(t) = Cx(t) + Du(t) = CT -1 x(t) + Du(t)Let Ā = TAT -1 , B = TB, C = CT -1 , then, after this transformation, we obtain the new state-space model . x = Ā x + Bu y = C x + Du Note that the inputs and outputs were not affected by this transformation; only the internal state vector and matrices changed. The transfer function corresponding to this model is given byH(s) = C(sI -Ā) -1 B + D = CT -1 (sI -TAT -1 ) -1 TB + D = CT -1 (sTT -1 -TAT -1 ) -1 TB + D = CT -1 T(sI -A) -1 T -1 TB + D = C(sI -A) -1 B + D = H(s)

Definition 6 .

 6 10 (A, B, C, D) and ( Ā, B, C, D) are algebraically equivalent if they are of the same dimension and they are related by a similarity (coordinate) transformation, i.e., there exists an invertible matrix T such that Ā = T -1 AT, B = T -1 B, C = CT, D = D.

Proposition 6 .

 6 11 Controllability, observability and stability are invariant under similarity transformations PROOF. Direct computations establish the following relations between the Gramians and State Transition M = T -1 M(T t ) -1 , N = T t NT, e Āt = T -1 e Āt T, matrices of two similar realizations.

.x-

  Ax = Bu y -Cx = Du in terms of the differentiation operator ∆ we obtain ∆ -

Definition 6 . 1 .

 61 18 (A, B, C, D) is reducible (non-minimal) if there exists a zero-state equivalent representation of smaller state-space dimension; otherwise (A, B, C, D) is irreducible or minimal. Some properties of minimal realizations are described by the following theorem. They make use of the so-called Hankel matrixH A,B,C = O(A, C)C(A, B) (A, B,C, D) is minimal iff it is controllable and observable. 2. (A, B, C, D) is minimal iff its Hankel matrix H A,B,C has full rank.

3 .

 3 Suppose (A, B, C, D) is minimal and let ( Ā, B, C, D) be a input-output equivalent realization. Then ( Ā, B, C, D) is minimal iff it is algebraically equivalent to (A, B, C, D). In such a case, the similarity transformation relating the two can be computed asT = C(A, B)C t ( Ā, B)(C( Ā, B)C t ( Ā, B)) -1 , T -1 = (O t ( Ā, B)O( Ā, B)) -1 O t (A, B)O( Ā, B).PROOF. 1. and 2. Let n be the dimensional of a minimal realization, n ≤ n. We note that rank(H A,B,C ) ≤ n, because if If n = n, follows from Sylvester inequality because rank(C(A, B)) ≤ n and rank(O(A, C)) ≤ n that implies rank(H A,B,C ) ≤ n = n. If n < n, let ( Ā, B, C, D) be an ndimensional realization. Since Ā is n × n, by Cayley-Hamilton theorem Ā n, ..., Ān-1 are linear combinations of I, Ā, ..., Ā n-1 . Hence rank(H A,B,C ) ≤ n. Suppose that (A, B, C, D) is controllable then rank(C(A, B)) = n and if (A, B, C, D) is observable rank(O(A, C)) = n then by Syvester inequality n + nn ≤ rank(O(A, C)C(A, B)) = rank(H A,B,C ) ≤ n that implies rank(H A,B,C ) = n. On other hand, n ≤ rank(H A,B,C ) ≤ n ≤ n that implies n = n and the realization is minimal.

c 1 a 2 b 2

 2

Proposition 6 .Proposition 6 . 24

 6624 [START_REF] Clarke | Functional Analysis, Calculus of Variations and Optimal Control[END_REF] Let T be an invertible transformation, then one has T -1 AT = Assume that the matrix T is irreducible and the products b i c i are positives for i = 1, ..., n. Then, the matrix T is similar to a symmetric tridiagonal matrix. Therefore, its eigenvalues are real. PROOF. Consider the matrix D = diag[d 1 , ..., d n ] where d 1 = 1 and d 2 i = d i-1 c i b i for i = 2, ..., n. Then it is readily verified that

Remark 1 .b 1 a 2 d 2 b 2 d 2 b 1 a 2 d 2 b 2 d 2

 122 In the similar way that in the last proposition, under the same assumptions, is possible define a diagonal positive matrix D = diag[d 1 , ..., d n ] where d 1 = 1 andd i = d i-1 b i c i (i = 2, ..., n), such that DT = H (equivalently, T = D -1 H)where H is a symmetric tridiagonal matrix determined from T as scaling[START_REF] Parlett | Reduction to tridiagonal form and minimal realizations[END_REF]) If Q -1 BQ = T is tridiagonal matrix and D is diagonal and invertible, then (QD) -1 B(QD) = D -1 TD is another tridiagonal matrix similar to B. T and D -1 TD are equivalent for theoretical purposes. Note that T i,i and T i+1,i T i,i+1 , i = 1, ..., n are invariant under diagonal scaling.

Figure 6 . 2 .

 62 Figure 6.2. An n-Compartment Catenary System (source [2]).

Figure 6 . 3 .

 63 Figure 6.3. An n-Compartment Mammilary System (source [2]).

  Solutions of the differential equation (2.9) for microalgae C. reinhardtii with parameter values µ = 2.34 [d -1 ], ρ = 0.34 [d -1 ], c = 0.253 [g.L -1]. considering regular intervals of 12 hours (0.5 [d]) and different initial conditions. A stable set exist in the interior of the interval [2.35, 2.8] over biomass concentration axis. . . . . . . . . . 2.4 Surface and level curves of the net rate of production (mean biomass volumetric productivity) for the optimization problem (2.6) for C. reinhardtii at the parameter values µ = 2.34 [d -1 ], ρ = 0.34 [d -1 ], c = 0.253 [g.L -1 ]. . . . . . . . . . . . . . . . . . . 2.5 Surface and level curves of the net rate of production for the optimization problem (2.6) for C. reinhardtii with the same parameter values in summer period. . . . . . . . 2.6 Surface and level curves of the net rate of production for the optimization problem (2.6) for C. reinhardtii with the same parameter values in regular time intervals. . . . xvi 2.7 Surface and level curves of the net rate of production for the optimization problem (2.6) for C. reinhardtii with the same parameter values in winter period. . . . . . . . . 2.8 Biomass trajectories associated to the optimal values of initial concentration x * 0 and terminal time T * in different light environments, with turnaround times t a = 1 (left) and t a = 2 (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Monod growth function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Phase plane of the differential equation (2.7) for microalgae C. reinhardtii with parameter values µ = 2.34 [d -1 ], ρ = 0.34 [d -1 ], c = 0.253 [g.L -1 ]. and different initial conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3

  .5. . . . . . . 3.3 Graph of the non-competitive inhibition function φ 2 (.) for the parameter values in table 3.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Growth functions with parameter values of table

  and parameters and functions involved are indicates in table 2.1.

Table 2 . 1 .

 21 Parameters and functions used in the general model.

	Parameter/function	Meanings

  2.7 Stability results a) If μ ≥ cρ then the differential equation (2.9) has stable solutions. b) If μ < cρ then the origin is a stable stationary point of the differential equation (2.9). PROOF. a) Is necessary to prove that a stable set (not necessarily stationary set) exist. The hypothesis μ ≥ cρ implies that dx dt > 0, ∀t ∈ [t 2k , t 2k+1 ], ie x(t) is increasing on this interval. Similarly, dx dt < 0, ∀t ∈ [t 2k+1 , t 2k+2 ], ie x(t) is decreasing on this interval for k = 0, 1, ..., this together with the continuity of x(t) implies that x(t 2k+1 ) is the maximum value of the function x in the interval [t 2k , t 2k+2 ]. If we define the sequence {x(t 2k+1 )} k∈N this sequence is bounded, then Bolzano-Weierstrass theo- rem says that there exists a convergent subsequence { x(t 2k+1 )} k∈N . Let U = lim k→∞ x(t 2k+1 ). Using a similar argument, we have that x(t 2k+2 ) is the minimum value of the function x in the interval [t 2k+1 , t 2k+3 ] and we can define the sequence {x(t 2k+2 )} k∈N that is bounded and there exists a convergent subsequence { x(t 2k+2 )} k∈N . Let L = lim k→∞ x(t 2k+2 ). It is clear from the definition of the sequences that U > L.

Table 2 .

 2 2. optimal values in different light environments for C. reinhardtii considering constant light.

Table 2 .

 2 2 sumarizes the results above for different light environments and different turnaround time values.

	constant light	0.2	2.7	3.674	0.9389
	Summer period [14-10 h]	0.3	1.6	1.91	0.6188
	Regular intervals [12 h]	0.4	1.5	1.71	0.5340
	Winter period [10-14 h]	0.8	0.4	1.4083	0.4345
	turnaround time	t a = 2
	Light environment	x *
	constant light	0.1	3.5	4.2475	0.7541
	Summer period [14-10 h]	0.2	2.6	2.4017	0.4786
	Regular intervals [12 h]	0.2	2.5	1.992	0.3981
	Winter period [10-14 h]	0.2	2.4	1.5743	0.3123

0 [gL -1 ] T * [d] x(T * ) [gL -1 ] J * [gL -1 d -1 ]

Table 2 . 3 .

 23 Parameter estimations for different incidental light settings for C. reinhardtii.

	I 0	µ	c
	1000 2.699 0.2736
	900	2.5918 0.2673
	800	2.4726 0.2604
	700	2.3395 0.2528
	600	2.1887 0.2444
	500	2.0147 0.2349
	400	1.8091 0.2239
	300	1.5576 0.2109
	200	1.2237 0.1948
	100	0.7778 0.1736

Table 2 . 4 .

 24 Optimal values in different incidental light settings for C. reinhardtii considering constant light.

	I 0	x * 0	T *	x (T * ) J *
	1000 0.2119 2.7270 4.3397 1.1075
	900	0.2062 2.7336 4.1663 1.0607
	800	0.2002 2.7424 3.9753 1.0087
	700	0.1936 2.7533 3.7620 0.9507
	600	0.1862 2.7674 3.5202 0.8849
	500	0.1777 2.7866 3.2411 0.8090
	400	0.1676 2.8142 2.9110 0.7193
	300	0.1554 2.8579 2.5068 0.6095
	200	0.1391 2.9442 1.9687 0.4639
	100	0.1148 3.1616 1.2454 0.2717
	this level.

Table 2 .

 2

	Parameter Value	units
	µ L	0.13 0.03	h -1 m
	K l Ī0	100 700	µmol. m -2 .s -1 µmol. m -2 .s -1

5.

Parameter values for C. reinhardtii obtained from

[START_REF] Janssen | Efficiency of light utilization of chlamydomonas reinhardtii under mediumduration light/dark cycles[END_REF]

.

Table 2 . 6 .

 26 Linear regression f (x) = axL

	µ [h -1 ]	x [g.L -1 ] I abs	I out	axL
	0.018	0.78	88	12	2.1203
	0.019	0.84	87	13	2.0402
	0.031	0.41	80	20	1.6094
	0.034	0.39	73	27	1.3093
	0.052	0.21	51	49	0.7134
	0.061	0.11	36	64	0.4463
	0.064	0.1	31	69	0.3711
	R 2				0.971
	Slope				2.33
	Zero-intercept a (estimated)	[m 2 . g -1 ]			2.81 0.0933

Table 2 . 7 .

 27 Parameters estimated for C. reinhardtii in this model.

	parameter	value	units
	μ ρ c	2.34 0.34 0.253	d -1 d -1 g.L -1

Table 3 . 1 .

 31 Variables used in the algal pond model

	Variable IC	Meaning
	x 1	1	Nitrifying biomass [gCOD/m 3 ]
	x 2	10	Microalgal biomass in C [gCOD/m 3 ]
	s 1	10	Ammonium concentration [gN/m 3 ]
	s 2	1	Nitrate concentration [gN/m 3 ]

  The equality µ 1 (s 1 ) = D give us an explicit expression for s 1 in this case, i.e., s * 1 = λ n . It exists if and only if D < µ 1 (s in ) (see Lemma 3.3, part 2).

Table 3 . 3 .

 33 Local stability conditions of non-coexistence equilibrium points for break-even concentrations

	Equilibrium point	Locally asymptotically stable if

Table 3 . 4 .

 34 Equilibrium existence and local stability conditions for break-even concentrations in system(3.18).

	Equilibria	Existence	Local stability
	Ē0	Always

Table 3 . 5 .

 35 Parameter values for the system(3.6). Graph of the non-competitive inhibition function φ 2 (.) for the parameter values in table 3.5.

	parameter	value	unit
	k 1	4.2	gN/gCOD
	k 2	0.08	gN/gCOD
	K n	0.7	gN/m 3
	K s	0.05	gN/m 3
	K p	0.1	gN/m 3
	K I	140	W/m 2
	a 0	0.2	1/m
	a 1	0.088	m 2 /gCOD
	μn	1.12	1/d
	μs	2.5	1/d
	μp	2.5	1/d

  • • • 5), as depicted on Figure 4.4 with the following diffusive exchange rates d 12 = 1 , d 13 = 2 , d 34 = 1 , d 35 = 3 , d 45 = 1 .

	d =1 12
	000000000 000000000 000000000 000000000 000000000 111111111 111111111 111111111 111111111 111111111
	d =2 13
		d =1 34
	d =3 35	d =1 45

  ds 2 +d(a 21 +a 23 +a 32 )s+da 21 a 32 s 3 +(a 12 +a 21 +a 23 +a 32 +d)s 2 +(a 12 a 23 +a 12 a 32 +a 21 a 32 +d(a 21 +a 23 +a 32 ))s+da 21 a 32

  Systems are input-output equivalent if and only if their transfer functions H C 4.4 and H T 4.6 are identical, that amounts exactly to have the four equalities, whatever is d. 32 = β 12 β 31 + β 13 β 21 β 12 + β 13 and replacing α 23 + α 32 by this last expression in (4.7) gives α 21 = β 21 + β 31 -β 12 β 31 + β 13 β 21 β 12 + β 13 13 (β 21β 31 ) 2 (β 12 + β 13 )(β 12 β 21 + β 13 β 31 ) So, we conclude that for any MRMT configuration with non-null parameters, there exists exactly one MINC equivalent that is input-output equivalent. Conversely, let S = β 21 + β 31 and P = β 21 + β 31 . β 21 and β 31 are solutions of X 2 -SX + P = 0. From (4.7) and (4.8), one has S = α 21 + α 23 + α 32 and P = α 21 α 32 . Then one has ∆ = S 2 -4P = (α 21 + α 23 + α 32 ) 2 -4α 21 α 32 = (α 21α 32 ) 2 + (α 23 ) 2 + 2α 21 α 23 + 2α 23 α 32 ≥ 0 So, there always exits two real non-negative roots, and accordingly to Remark 1, there are two symmetric solutions for β 21 and β 31 . Assume without loss of generality that β 21 ≥ β 31 : β 21 = α 21 + α 23 + α 32 + (α 21 + α 23 + α 32 ) 2 -4α 21 α 32 2 β 31 = α 21 + α 23 + α 32 -(α 21 + α 23 + α 32 ) 2 -4α 21 α 32 2 When α 23 > 0, note that β 21 and β 31 are distinct. Then, there exits a unique pair (β 12 , β 13 ): 12 α 21α 23α 32 + (α 21 + α 23 + α 32 ) 2 -4α 21 α 32 2 (α 21 + α 23 + α 32 ) 2 -4α 21 α 32 β 13 = α 12 -α 21 + α 23 + α 32 + (α 21 + α 23 + α 32 ) 2 -4α 21 α 32 2 (α 21 + α 23 + α 32 ) 2 -4α 21 α 32 Finally, note that the property (α 21 + α 23 + α 32 ) 2 -4α 21 α 32 -(α 21α 23α 32 ) 2 = 4α 21 α 23 ≥ 0 implies that the expressions of β 12 and β 13 given above are non negative.

	α 21 + α 23 + α 32 = β 21 + β 31	(4.7)
	α 21 α 32 = β 21 β 31	(4.8)
	α 12 = β 12 + β 13	(4.9)
	α 12 (α 23 + α 32 ) = β 12 β 31 + β 13 β 21	(4.10)
	Equation (4.9) gives immediately: Replacing α 12 in (4.10) gives α 23 + α or equivalently α 12 = β 12 + β 13 α 21 = β 12 β 21 + β 13 β 31 β 12 + β 13 Then, condition (4.8) provides α 32 = β 21 β 31 (β 12 + β 13 ) β 12 β 21 + β 13 β 31 Finally, (4.10) gives α 23 = β 12 β 31 + β 13 β 21 β 12 + β 13 -β 21 β 31 (β 12 + β 13 ) β 12 β 21 + β 13 β 31 or equivalently α 23 = 1 1 β 31 β 21 β 12 β 13 = α 12 1 α 23 + α 32 β 12 β 13 = α 12 β 21 -β 31 β 21 -1 -β 31 1 1 α 23 + α 32 that are β 12 β From (4.9), (4.10), one can write β 12 = α

This chapter is based on the paper H. Ramírez C., A. Rojas-Palma and D. Jeison, Productivity optimization of microalgae cultivation in a batch photobioreactor process, submitted to Mathematical Methods in the Applied Science, Wiley (2016).

Remark Note that the last result does not fully characterize the nature of the coexistence equilibrium but at least it avoids the possibility of existence of Hopf bifurcations and other local bifurcations for such equilibrium point for the different parameter values of the model.

This chapter is based on the paper A. Rapaport, A. Rojas-Palma, J. R. de Dreuzy and H. Ramírez C., Equivalence of finite dimensional input-output models of solute transport and diffusion in geosciences, submitted to IEEE Transactions on Automatic Control, IEEE Xplore (2016).

• (µ 1 (s 1 ) -D)x 1 = 0 ⇒ x 1 = 0 or µ 1 (s 1 ) = D,

• (µ 2 (s 1 , x 2 ) + µ 3 (s 1 , s 2 , x 2 ) -D)x 2 = 0 ⇒ x 2 = 0 or µ 2 (s 1 , x 2 ) + µ 3 (s 1 , s 2 , x 2 ) = D.

1. If x 1 = 0 and x 2 = 0, then ṡ1 = 0 ⇒ D(s in -s 1 ) = 0 ⇒ s 1 = s in , ṡ2 = 0 ⇒ -Ds 2 = 0 ⇒ s 2 = 0. Therefore (x 1 , x 2 , s 1 , s 2 ) = (0, 0, s in , 0) is an equilibrium point.

2. If µ 1 (s 1 ) = D and x 2 = 0, let λ n = µ -1 1 (D) the unique solution of the first equality, then

Therefore (x 1 , x 2 , s 1 , s 2 ) = s in -λ n k 1 , 0, λ n , s inλ n is an equilibrium point. In particular, we can see that µ 1 (s 1 ) is an increasing function and D is constant, then to ensure the intersection between these functions it is necessary that D < µ 1 (s in ), the last inequality is a necessary and sufficient condition to guarantee the existence of this equilibria.

3. If x 1 = 0 and µ 2 (s 1 , x 2 ) + µ 3 (s 1 , s 2 , x 2 ) = D, suppose that s 2 = 0, the last equality is reduced to µ 2 (s 1 , x 2 ) = D and ṡ1 = 0 ⇒ D(s in -s 1 ) -k 2 µ 2 (s 1 , x 2 )x 2 = 0 ⇒ x 2 = s in -s 1 k 2 , if we denote λ s the unique positive solution of the equation

then (x 1 , x 2 , s 1 , s 2 ) = 0, s in -λ s k 2 , λ s , 0 is an equilibrium point. Similarly to the previous case, µ 2 (s 1 , s in -s 1 k 2 ) is an increasing function and D is constant, then to ensure the intersection between these functions it is necessary and sufficient that D < µ 2 (s in , 0). From the analysis of the above cases, the statement is proved.

Model reduction

The dynamics described by the system (3.6) can be studied from a reduced model based on some simple assumptions. Let us define:

whose dynamics writes: ż = -Dz.

Then the system (3.6) is equivalent to the following

(3.17)

Where ż = -Dz, we have z(t) = z 0 e -Dt with z 0 = z(0) ≤ 1 by (3.16) and obviously lim t→+∞ z(t) = 0. Hence, using the theory of asymptotically autonomous systems [START_REF] Smith | The Theory of the Chemostat[END_REF][START_REF] Strauss | On asymptotically autonomous differential equations[END_REF][START_REF] Thieme | Convergence results and a poincaré-bendixson trichotomy for asymptotically autonomous differential equations[END_REF] one studies the limiting system (we can consider the system (3.6) restricted to the invariant hyperplane z = 0) (3.18) where, from (3.16) µ 3 (s 1 , x 2 ) = µ 3 (s 1 , s in -k 2 x 2 -s 1 , x 2 ),

This system is defined in the region

and all the parameters are positive.

There are simple hypotheses to be checked before one can conclude that the dynamics of the original system (3.6) and that of the asymptotic limiting equations (3.18) have the same asymptotique behavior [START_REF] Hsu | A survey of mathematical models of competition with an inhibitor[END_REF], see Appendix F in [START_REF] Smith | The Theory of the Chemostat[END_REF]. The important hypothesis is the lack of a cyclic connection for orbits, and we comment on this since it is an important hypothesis for uniform persistence. We describe only the case of equilibrium points, the stability of them prevents it from being a part of any chain of equilibria (we see this below).

Remark

The function µ 3 (.), defined in the region

is positive, smooth, bounded, decreasing with respect to s 1 and x 2 and vanishes over the line

we obtain a bound for the norm of ∇V • f 2 (x 1 , x 2 , s 1 )

,

where

From the last inequalities

. Therefore Theorem 4.18 in [START_REF] Khalil | Nonlinear Systems: Pearson New International Edition. Always learning[END_REF] guarantee that for all initial conditions (x 1 (t 0 ), x 2 (t 0 ), s 1 (t 0 )) ∈ R 2 ⊂ Ω 2 , there exist constants k, γ, b, T, such that ∀ t 0 ≤ t ≤ t 0 + T,

and

. This completes the proof.

Remark

We have that the system (3.6) is dissipative (see remark of Lemma 3.2), all the finite equilibrium points are hyperbolic and the system (3.17) does not possess a cycle of rest points.

) is a solution of the system (3.17) then lim

where system (3.17) is equivalent to the system (3.6), all solutions of this system with initial condition in R ⊂ Ω tend asymptotically to the positive coexistence equilibrium point E c in a finite period of time, and then remain bounded and near equilibrium, following the result of Theorem 3.10.

ordinary differential equations:

where the parameters d ij = d ji (i = j) denote the diffusive exchange rates of solute between reservoirs i and j.

For the sake of simplicity, we shall assume

which is always possible by a change of the time scale of the dynamics.

In the following we adopt an input-output setting in matrix form:

where X denotes the vector of the concentrations and the matrix A satisfies the following properties.

Assumptions 4.1 There exist matrices V and M such that

where B is defined in (6.8), V is a positive diagonal matrix and M is a symmetric matrix that fulfills i. M is irreducible (i.e. the graph with nodes P i , and edges -→ P i P j when M ij = 0 is strongly connected)

ii. M ii > 0 for any i

Let us decompose the matrix M as follows

where L is a row vector of length n -1. Then equality MX = -V 1 X 1 B amounts to write

M being invertible (Lemma 4.4), one can write X = -M-1 L ′ X 1 and thus X 1 has to fulfill

From Assumptions 4.1, one has M✶ = 0 which gives

We conclude that one should have X 1 = 0 and then X = 0, that is X = 0. The matrix A is thus invertible.

Finally the system admits a unique equilibrium X ⋆ = -A -1 Bu for any constant control u. As Assumptions 4.1 imply the equality A✶ = -B, we deduce that the equilibrium is given by X ⋆ = ✶u.

Lemma 4.6 Under Assumptions 4.1, the sub-matrix à is diagonalizable with real negative eigenvalues.

PROOF. Note first that the matrix à can be written as à = -Ṽ-1 M. The matrix Ṽ being diagonal with positive diagonal terms, one can consider its square root Ṽ1/2 , defined as a diagonal matrix with Ṽi terms on the diagonal, and its inverse Ṽ-1/2 . Then, one has

which is symmetric. So à is similar to a symmetric matrix, and thus diagonalizable. Let λ be an eigenvalue of Ã. There exists an eigenvector X = 0 such that

As M is definite positive (Lemma 4.4) as well as Ṽ, we conclude that λ has to be negative.

About controllability and observability

The controllability and observabillity definitions for single-input single-output systems (A, B, C) of dimension n are given in part B of Annexes. We consider now the following definition Definition 4.7 To a given triplet (A, B, C), we associate the linear operator F A,B,C : L

for the initial condition X(0) = 0. We say that a triplet (A, B, C) is a minimal representation if among all the triplets

We recall a well known result of the literature on linear input-output systems [59, 2.4.6]. The particular structures of the matrices A, B and C that we consider allow to show the following property. Lemma 4.9 Under Assumptions 4.1, (A, B) controllable is equivalent to (A, C) observable.

PROOF. Note first that one has

Then one can write

The Multi-Rate Mass Transfer and Multiple INteracting Continua configurations

We consider two particular structures of networks whose (A, B, C) representations fulfill Assumption 4.1.

where the immobile parts are connected in series, of length n -1, one of them being connected to the mobile zone (see Fig. 4.3). In the context of general compartmental models, the flow digraph asociated to matrix A is an n-compartment catenary system [2].

In the following, we give properties on eigenvalues for MRMT matrices only, because it is easier to be proved for this particular structure. In the next section, we shall show that MINC and MRMT structures are indeed algebraically equivalent, and as a consequence eigenvalues of MINC matrices fulfill the same properties. Lemma 4.12 A MRMT matrix is Hurwitz (i.e. all the real parts of its eigenvalues are negative).

PROOF. Take a number

Then the matrix γI + A is an irreducible non-negative matrix. From Perron-Frobenius Theorem (see [49, 8.4.4]), r = ρ(γI + A) is a single eigenvalue of γI + A and there exists a positive eigenvector associated to this eigenvalue. That amounts to claim that there exists a positive eigenvector X of the matrix A for a single (real) eigenvalue λ = rγ, and furthermore that any other eigenvalue

From the particular structure of MRMT matrix, such a vector X has to fulfill the equalities

The vector X being positive, we deduce that λ is negative.

Theorem 4.17 Any minimal representation (A, B, C) that fulfills Assumptions 4.1 is algebraically equivalent to a MRMT structure.

Remark For a positive linear system (A, B), let A + 0 (A, B) be the attainability set from the 0-state with non-negative controls. The system being positive, one has A + 0 (A, B) ⊂ R n + and for any state X ∈ A + 0 (A, B), the state Z = R -1 X for the equivalent MRMT structure is also non-negative, but for a state X ∈ R n + \ A + 0 (A, B), the equivalent state Z = R -1 X is not necessarily non-negative (as the coefficients of the matrix R -1 are not necessarily non-negative). Consequently, one can have an equivalent input-output representation in MRMT form but with negative concentrations.

Getting back to the original system (4.1) we note that it is a classical linear time-invariant (LTI) single input-single output (SISO) model. The transfer matrix associated to this system is defined by

this matrix function relates the input and output of the system and some results can be determined from this expression (stability, poles, minimal realization for instance). If the system (A, B, C) fulfills the assumption 4.1 and is controllable, then it is also observable (by Lemma 4.9) and the realization (A, B, C) is minimal. Also

then the systems (A, B, C) and (A MRMT , B, C) are transfer equivalents (See appendix B in Annexes). Therefore, from Theorem 4.17, the equality above and the minimality of the realizations we have the following result Corollary 4. [START_REF] Chen | Linear System Theory and Design[END_REF] The minimal realizations (A, B, C) and (A MRMT , B, C) are input-output equivalents.

Equivalence with MINC structure

Take a matrix A that fulfills Assumption 4.1 and such that pair (A, B) is controllable. As we have already shown that such representation (A, B, C) is minimal and equivalent to a MRMT configuration, we can assume without any loss of generality that the matrix A has the structure

where ∆ is a square diagonal matrix (of size n -1) with distinct negative eigenvalues. We denote by V the diagonal matrix of the volumes associated to the matrix A with V 1 = 1, as given by Lemma 4.2 We shall consider a tridiagonalization of this matrix. For this purpose, we recall the Lanczos algorithm.

with the corresponding transfer function TMINC (z) = 6 + 8.6z + z 2 6 + 35.4z + 12.6z 2 + z 3 .

The Nyquist plots of the transfer functions T, TMRMT and TMINC are reported on Figure 4.5, showing the quality of the approximation with only three compartments derived from the MRMT or MINC representations. There exist many reduction methods in the literature, but a reduction through MRMT or MINC has the advantage to obtain easily reduced models with a physical meaning. Remark For a positive linear system (A, B), let A + 0 (A, B) be the attainability set from the 0-state with non-negative controls. The system being positive, one has A + 0 (A, B) ⊂ R n + and for any state X ∈ A + 0 (A, B), the state Z = R -1 X for the equivalent MRMT or MINC structure is also nonnegative, but for a state X ∈ R n + \ A + 0 (A, B), the equivalent state Z = R -1 X is not necessarily nonnegative (as the coefficients of the matrix R -1 are not necessarily non-negative). Consequently, one can have an equivalent input-output representation in MRMT form but with negative concentrations for such states of the system. Theorems 4.17 and 4.26 show that whatever the network structure is, it is always possible to represent its input-output map with either a MRMT star structure or a MINC series structure. Nevertheless, these two results require that the original network is of minimal representation, or equivalently that the system (4.1) is controllable (or observable, see Proposition 4.9). To illustrate the necessity of the controllability assumption, we present an example of a structure of four reservoirs of volumes (see Fig. 4.6)

Im

with the diffusive exchange rate coefficients:

equivalent matrix A MRMT in a MRMT form such that that fulfills the assumption 4.1, then there exist a diagonal positive matrix V and a symmetric matrix M such that

The next result determine that the trace of matrices V and V are invariant by transformation.

Proposition 4.27 Under the hyphotesis of Theorem 4.17, the similarity transformation R preserve the trace of the matrices V and V, ie, tr(V) = tr( V).

PROOF. In first place, where

From the last equality, we have

with this tr(V) = -B t A -1 ✶ n , and using the same arguments we obtain tr( V) = -B t A -1 MRMT ✶ n , and

Therefore tr(V) = tr( V) and the transformation R preserve the traces of V and V.

Remark In the proposition above, using the same argument but for a tridiagonal matrix A MI NC that fulfills the Assumption 4.1, then exist a diagonal positive matrix V and a symmetric tridiagonal matrix M such that

So, the result above is also valid for a MINC instead a MRMT structure.

Chapter 6

Annexes A Linear Algebra results Proposition 6.1 Let A ∈ M n be positive semi-definite and X ∈ R n . Then X t AX = 0 iff Ax = 0.

PROOF. Suppose that X = 0 and X t AX = 0. Consider the polynomial p(t)

The hyphotesis assure that p(t) ≥ 0 for all real t. However, if AX 2 = 0 then for sufficiently large negative values of t we would have p(t) < 0. We conclude that AX 2 = 0, so AX = 0. Corollary 6.2 A positive semi-definite matrix is positive definite iff it is nonsingular (invertible).

PROOF. Suppose that A is positive semi-definite. Now, A singular implies AX = 0 for some X ∈ R n no null, but then X t AX = 0 for some X ∈ R n and A is not positive definite.

Proposition 6.3

Let A ∈ M n Hermitian (symmetric in case of K = R) and C ∈ M n,m . Suppose that A is positive definite, Then rank(C t AC) = rank(C) and C t AC is positive definite iff rank(C) = m.

PROOF. For proof and details see [START_REF] Horn | Matrix Analysis. Matrix Analysis[END_REF], page 431, observation 7.1.8. PROOF. To prove positive (semi-) definiteness of G, we need to examine the associated quadratic form q(X) = X t GX = X t A t AX = (AX) t AX = AX 2 ≥ 0, for all X ∈ R n . Moreover, it equals 0 if and only if AX = 0, and so if A has trivial kernel, this requires X = 0, and hence q(X) = 0 if and only if X = 0. Thus, in this case, q(X) and G are positive definite.

Proposition 6.6 1. The determinant of a Gram matrix is equal to zero iff any of its principal minors is zero.

A real matrix and its Gram matrix have the same rank.

PROOF.

1. See for instance [START_REF] Eves | Elementary Matrix Theory[END_REF].

2. Let A real matrix, G = A t A their Gram matrix asociated and X vector such that GX = 0. Then X t A t AX = 0 and this implies AX 2 = 0. We therefore have AX = 0. Conversely, if X satisfies the relation AX = 0, then clearly A t AX = 0. It folllows that the homogeneous systems AX = 0 and GX = 0 are equivalent. Hence, by dimensionality theorem for homogeneous systems

where n is the number of columns of A. The theorems is therefore proved.

Theorem 6.7 Perron-Frobenius

Let A be an n × n matrix with nonnegative real entries and irreducible. Then we have the following 1. A has a nonnegative real eigenvalue. The largest such eigenvalue, λ A , dominates the absolute values of all other eigenvalues of A. The domination is strict if the entries of A are strictly positive.

2. If A has strictly positive entries, then λ A is a simple positive eigenvalue, and the corresponding eigenvector can be normalized to have strictly positive entries.

3. If A has an eigenvector v with strictly positive entries, then the corresponding eigenvalue λ v is λ A .

PROOF. See [START_REF] Horn | Matrix Analysis. Matrix Analysis[END_REF] for details.

B Realization Theory Fundamentals

The state-space model of a continuous-time linear system is given by .

x = Ax + Bu y = Cx + Du (6.1)

114

The vector x is called the state vector of the system. We will denote the number of states in the system by n, so that x ∈ R n . The quantity n is often called the order or dimension of the system. In general, we might have multiple inputs u 1 , u 2 , ..., u m to the system. In this case, we can define an input vector u = [u 1 , u 2 , ..., u m ] t . In the same way, we might have multiple outputs y 1 , y 2 , ..., y p . In this case, we can define the output vector y = [y 1 , y 2 , ..., y p ] t . Note that each of these outputs represents a sensor measurement of some of the states of the system.

The system matrix A is an n × n matrix representing how the states of the system affect each other. The input matrix B is an n × m matrix representing how the inputs to the system affect the states. The output matrix C is a p × n matrix representing the portions of the states that are measured by the outputs. The feedthrough matrix D is a p × m matrix representing how the inputs affect the outputs directly (i.e., without going through the states first).

While the state-space model is a time-domain representation of a system, one can also convert them to the frequency domain by taking the Laplace transform in continuous-time. Specifically, if we take the Laplace transform of (6.1), we obtain:

Note that this includes the initial conditions of all the states. The first equation can be rearranged to solve for L {x} (s) as follows:

Substituting this into the equation for L {y} (s), we obtain

Definition 6. [START_REF] Ballyk | Global analysis of competition for perfectly substitutable resources with linear response[END_REF] The transfer function of the state-space model (6.1) with initial condition x(0) = 0 is

Note that H(s) is a p × m matrix, and thus it is a generalization of the transfer function for standard single-input single-output systems. In fact, it is a matrix where entry i, j is a transfer function describing how the jth input affects the ith output.

The basic problem in Realization Theory is to determine matrices A, B, C, D such that the state space system (A, B, C, D) of system (6.1) represents a given input-output map, specified by its impulse response or transfer function.

We recall the usual definitions of controllability and observability of single-input single-output systems (A, B, C) of dimension n (see for instance [START_REF] Kailath | Linear Systems. Information and System Sciences Series[END_REF]). Definition 6.9

• The controllability matrix associated to the pair (A, B) is given by

Definition 6.12 Two systems of the form (6.3) are called transfer equivalent (or zero-state equivalent) if and only if their transfer matrices are identical.

Obviously, transfer equivalence is an equivalence relation in the set of all systems of the form (6.3). Definition 6.13 Two systems of the form (6.3) are called input-output equivalent if they are satisfied by the same pairs (x, u), ie, the systems has to be at zero state and given a L 2 signal u, the output signals y should be identical.

Is clear from the last definition that two system input-output equivalent has the same inputoutput behaviour. Let us recall that a square polynomial matrix K(s) is called unimodular if detK(s) is constant and different from 0. Then the inverse of K(s) is also a polynomial unimodular matrix. It is known [START_REF] Wolovich | Linear Multivariable Systems[END_REF] that K(s) is unimodular if and only if it can be obtained from the identity matrix by finitely many elementary row operations over the ring of polynomials in s. The next results give us criteria to determine when two system are input-output and transfer equivalents. Their proofs are in [START_REF] Blomberg | Algebraic theory for multivariable linear systems[END_REF]. Proposition 6.14 Two input-output systems given by matrices [P 1 (s), Q 1 (s)] and [P 2 (s), Q 2 (s)] are transfer equivalent if and only if there are polynomials matrices M 1 (s), M 2 (s)

Proposition 6.15 Two input-output systems given by matrices [P 1 (s), Q 1 (s)] and [P 2 (s), Q 2 (s)] are inputoutput equivalent if and only if there is a unimodular matrix K(s) such that

Corollary 6. [START_REF] Carrera | On matrix diffusion: formulations, solution methods and qualitative effects[END_REF] If two systems are input-output equivalent then they are transfer equivalent. Now, consider the system (6.1) and .

x = Ā x + Bu y = C x + Du (

we call these systems (A, B, C, D) and (A, B, C, D) respectively. When will these two systems have the same transfer function? Theorem 6.17

Lemma 6. [START_REF] Clarke | Nonsmooth Analysis and Control Theory[END_REF] The unique solution of the system TX = B with T invertible and irreducible tridiagonal matrix has all its entries non zero.

PROOF. In first place, we note that the system TX = B with X = (x 1 , . . . , x i-1 , x i , x i+1 , . . . , x n ) t , can be write as

Suppose that in the last system x 1 = 0, then

wich clearly is a contradiction. Then necessarily x 1 = 0 and by the same argument X(2 : n) = 0 n-1 .

In fact if x n = 0 then x n-1 = 0 (T is tridiagonal) and a recursive argument implies that X(2 : n) = 0 n-1 which contradicts the last result, reason why x n = 0 (and in consequence x n-1 = 0). Suppose now that for some i ∈ {2, . . . , n -2}, x i = 0. Also, where T is no singular, the system TX = B has unique solution.

with x i = 0 and y i = 0 (i = 1 • n) if and only if

Corollary 6.23 Let H n the n × n Hankel matrix corresponding to the system (A, B, C) where n is the dimension of A and S be an invertible transformation, where H = S -1 AS is tridiagonal irreducible (ie, with non null entries in the first diagonal below and the first diagonal above the principal diagonal), S -1 B = x 1 B and CS = y 1 C with x 1 , y 1 = 0. There exists a transformation of (A, B, C) to a system (S -1 AS, S -1 B, CS) if and only if all the leading principal minors of H n are non-zero.

Note that in proposition above T -1 C(A, B) needs to be upper triangular. The next algorithm allows to construct a matrix L such that LC(A, B) be upper triangular with positive elements in the diagonal, then we can consider L = T -1 the transformation of the statement.

Algorithm 2 Upper triangular matrix with positive diagonal elements

Choose matrix A ∈ M n and I ∈ M n identity matrix. Consider Ā = [A|I] and E i (-1) elementary matrix. For i = 1, ..., n -1 If A(i, i) = 0, For j = i + 1, ..., n Ā(j, :) = Ā(j, :) -Ā(i, :) * Ā(j, i)/ Ā(i, i),

End For End If

End For P = Ā(1 : n, n + 1 : 2n), where the amounts of tracer, x i (t), with i = 1, 2, ..., n, are referred to as the state variables of the system, the matrix A is called the system matrix, and the input vector b(t) is known as the forcing function.

Many properties of model (6.12) depend only on the disposition of zero and nonzero elements in the compartmental matrix A (matrix singularity of A turns out to be such a property). If a ij = 0 with i = j, then there is no flow of material from compartment j directly to compartment i. If a ij = 0 (i = j), then material is passing directly from the j th to the i th compartment. A useful representation of structure of a compartmental system is the connectivity diagram for the system which shows the nonzero transfers of material on a directed graph. The directed graph consists of a set of nodes (the compartments) together with a set of directed edges (the flows) connecting certain of these nodes. A compartmental system, complete with inputs and exits, and its associated connectivity diagram are given in 6.1.

General compartmental model

We consider now a network of n interconected compartments, and assume that a same single species is presented in each vessel at initial time. The network is fed with a single limitant ressource that allows th growth of the bacterial species.

The nodes of the interconnection graph is composed of the tanks of volume denoted by V i > 0 (i = 1, ..., n), and the arcs represent mass transfers between the tanks. We consider four types of arcs:

• Flux of flow rate Q ij ≥ 0 from node i to node j with j = i,

• Fick diffussion of parameter d ij = d ji ≥ 0 between nodes i and j with j = i

and assume that the Kirchoff law or mass conservation is satisfied at each nodes, wich amouts to write thee equalities

We denote by I and O the sets of input and output nodes

We shall represent by S and X respectively the vectors of concentrtions of the limiting ressource and biomass in the set of reactors. The time evolution of these vectors is modelled by the following dynamical system Ṡ = -1 y R(S, X) + MS + DS in Ẋ = R(S, X) + MX + DX in (6.15) where y > 0 and R(.) stand respectively for the convertion yield of the species, and the vector of kinetics that occurs in each vessel. Classicaly, we assume the kinetics to be linear w.r.t. to the biomass concentration and characterized by a specific growth rate µ(.) such that R i (S, X) = µ(S i )X i , ∀ i = 1, ..., n. (6.16)

Remark Under this last assumption, one can check that the parameter y can be chosen equal to one without any loss of generality, simply replacing X i /y by X i in equation (6.15).

We shall assume that the function µ(.) is smooth (at least C 2 ) with µ(0) = 0 and µ(S) > 0 for S > 0.

The matrices M and D represent respectively the mass transfers and inputs and are defined as follows: An matrix A such that -A is compartmental and diagonal dominant are called M-matrices.

Lemma 6.27

The matrix M is compartmental PROOF. The first two properties are fullfilled from (6.17) and for any i, one has

and using equality (6.13) we obtain

Thus propertiy 3. is fulfilled.

In the following, we assume that there is no loop in the network that is not connected to the set O. Assumption 6. [START_REF] Donado | Multicomponent reactive transport in multicontinuum media[END_REF] The sets I and O are non-empty and for i / ∈ O there exists j ∈ O and a sequence {i 0 , ..., i k } such that i 0 = i, i k = j with Q i α i α+1 > 0. Proposition 6.29 Assume that assumption 6.28 is fulfilled. The positive domain D = R 2n + is invariant by the dynamics (6.15) and any solution in D is bounded. Furthermore, any solution of the system (6.15) fulfills lim t→∞ S(t) + X(t) = -M -1 D(S in + X in ) (6.21) PROOF. From system (6.15) and the property that the matrix M has non-negative off-diagonal terms, one deduces the following inequalites

and cnclude that the domain D is invariant by the dynamics (6.15). Now consider the vector Z = S + X and for convenience we posit Z in = S in + X in . One has straithforwardly Ż = MZ + DZ in . (6.22) We recall that for a compartmental matrix, if the graph is outflow-connected implies that the matriz M is non-singular and Hurwitz. We can then deduce the asymptotic property of solution of system (6. Furthermore, a non-singular comparmental matrix being diagonal dominant, -M is a M-matrix. The components of the vector DZ in being non-negative, the properties of the M-matrices implies then the vector Z has non-negative components. Finally, we deduce the boundedness of the solutions S and X.

More detailed analysis for linear and nonlinear compartmental systems can be founded in [START_REF] Jacquez | Compartmental analysis in biology and medicine : Kinetics of distribution of tracer-labeled materials[END_REF][START_REF] Jacquez | Qualitative theory of compartmental systems[END_REF].

List of publications and conferences

Articles (unpublished)