- Email: [email protected]
Qualitative Dynamics of a Class of Nonlinear Biological Systems
Copyright 10 IFAC System Structure and Control, Nantes, France, 1998
QUALITATIVE DYNAMICS OF A CLASS OF NONLINEAR BIOLOGICAL SYSTEMS Olivier BERNARD· Jean-Luc GOUZE ** • GESAME, UGL , Av. G. Lemaitre 4, B-1348 Louvain-La-Neuve, Belgium, [email protected] ** GOMORE Project, INRIA , BP 93, 06902 Sophia-Antipolis, France, [email protected]
Abstract: In this paper we propose a methodology to analyze the global qualitative dynamics of a nonlinear class of differential systems which encompass numerous biological models. We show that, independently of parameters values, the possible succession with respect to time of some qualitative events characterizing transients of state variables (extrema, position ... ) are strongly related to the signs of the Jacobian matrix. The method is illustrated with a set of models usually used to describe phytoplanctonic growth in a bioreactor. The corresponding transition graph is derived and compared with experimental data. This methodology applies well when some parameters are unknown , which is often the case in biological models. Copyright © 1998 IFAC
Resume : Nous proposons une methode d 'analyse de la dynamique qualitative globale d'une classe de systemes differentiels non-lineaires, decrivant de nombreux modeles biologiques. On montre que la succession temporelle de certains evenements dynamiques qualitatifs (passage par un maximum ou un equilibre) depend des signes de la matrice jacobienne. On illustre la methode sur un modele de croissance du phytoplancton dans un bioreacteur. On trace le graphe de transition, qui est compare aux donnees experimentales. Cette methode ne depend pas de la valeur precise des parametres. Keywords: qualitative analysis, nonlinear systems, validation, ecology, biotechnology.
1. INTRODUCTION
systems with monotonous interactions) which has numerous applications in the biological field, including gene regulation models, compartmental systems (Goodwin, 1965j Levine, 1985), cellular growth (Burmaster, 1979), and development of stage structured populations (CasweU, 1989).
In this paper we propose to study the link between the structure of an autonomous differential system (defined by the signs of the Jacobian matrix and its zeros) and its dynamical behavior, in terms of qualitative features for the state variables. The main idea is to describe independently of parameters the dynamics allowed by the structure of a model in order to validate (or invalidate) this structure after a data analysis. In particular, the precise value of some parameters may be unknown, as it is often the case in biological models. We focus more particularly on a class of autonomous differential systems (loop structured
For this class of models, a methodology has recently been proposed (Bemard and Gouze, 1995bj Bemard and Gouze, 1995a) to derive the set of possible qualitative behavior for the state variables. The transient behavior is described either by their possible successive tendencies (increasing or decreasing) (Bemard and Gouze, 1995b) , or by their possible successive position toward
763
O"J E {-I , 1},q E {1 , ... , 2n }}
equilibrium value (larger or lower) (Bemard and Gouze, 1995a). From this analysis it is possible to test experimentally if the extradiagonal signs of the Jacobian matrix are consistent with the data.
Definition 3. For x* E 0 the system with monotonous interactions (E) is diagonally x*-monotonous if for all j in {I , ... , n} , for all q in {I , ... , 211. }, the sign of the partial derivative 8Ii/8xj(x) is fixed on each domain:
In this paper we generalize this result by considering a more precise qualitative state describing simultaneously the tendency of the state and its position toward equilibrium. We consider thus the domain 0 as partitioned by the different regions delimited by the nullclines and the hyperplanes associated to an equilibrium x*. We give in section 2 the theorems that describe the global qualitative behavior of the system from structural considerations i. e. the way trajectories go from one of these regions to another. We present also a graphical version of these theorems, easier to use for validation purposes. Finally, in section 3 we study a class of systems modeling growth of phytoplankton in a continuous bioreactor (chemostat) . The qualitative transient behavior is presented, and some results are derived. It is shown that at the most one maximum, one minimum, one equilibrium crossing bottom-up and one top-down for each state variable is possible. The transition graph obtained is then compared to some experiments, and conclusions on the structure of the models usually used to describe phytoplanktonic growth are drawn.
The domains WO". (x*) are called the orthants of the x*-deviation space or x*-orthants. In the same way, we define the orthants of the velocity space or z-orthants:
Definition 4. For (O"q , O"P) E S~ let us define the following open set: 00"'0"1' (x*) ~f WO" . (x*) n ZO"1'
{x E WO".(x*) ; diag(O"P)I(x)
> O}
Note that some 00"00"1' (x* ) are empty. The remaining non empty 00"' 0"1' (x*) represents therefore the qualitative situations allowed by the model and will be called the possible regions. It is shown in (Bernard and Gouze, 1997) how to derive these possible regions from the signs of the J acobian matrix. We remove from the 00"' 0"1' (x* ) some manifolds of zero measure containing trajectories that may have undesiderable behaviors, and we denote the corresponding non empty domains by nO".O"1'(x*) . We will then show that the transition between these domains obey some rules determined by the extradiagonal terms of the J acobian matrix.
Due to the lack of space, the proofs and details are to be found elsewhere (Bernard and Gouze, 1995b; Bernard and Gouze, 1997) .
2. DEFINITIONS Notations. The notations x > 0 for x = \Xl , ... ,x n ) E Rn means that for all i, Xi> O. For x ERn, sign(x) is the vector with components sign(xi} (with values taken in {-l , O, l}).
Let 0 be an open convex domain of Rn and
=
I a
3. THEOREMS
Cl mapping from 0 onto Rn. We consider on 0 the autonomous differential system :
Theorem 1. (Transitions between regions). Consider a loop system (E) with monotonous interactions and an equilibrium point x*. Suppose that nO"n 0"1'1 (x*) and n".n".1' 2 (x*) are two strict neighbors (i.e. the manifold which separates these regions is of dimension (n-l)) . We denote tk,k+l the sign of the (k, k + 1) element of the Jacobian matrix. A trajectory can go from nO" ' 2".1'2 (x*) to n".n ".1'1 (x*) (resp n""l ".1'1 (x*) to nO"Q2 ".1'2 (x*)] if one of the two following conditions is verified:
(E){x=I(x) Definition 1. A system (E) has a loop structure if hex) = h(Xi , Xi+r) Vi E {1, ... ,n}. Definition 2. The system (E) has monotonous interactions on 0 if each partial derivative 8 Id 8x j (x) for i :f. j never cancels on O.
= =
• tk , k+lO":~l 0":1 (resp. -0":1]; then Xk admits a minimum [resp. a maximum] . • tk,k+l 0"~~1 0"~1 [resp _0"~1]; then Xk crosses its equilibrium xA; bottom-up Iresp top-down].
The signs of the elements define what we call the structure of the system. We will consider the set Sn containing 2n elements:
We say that n""lO"1'l (x*) Iresp accessible from n".'2".1'2 (x*) Iresp 764
n".'2".1'2 nO"n".1'l
(x*)] is (x*)] .
The following theorem describes the behavior of the trajectories of a differential system (~), loop structured and diagonally-monotonous, with monotonous interactions.
4. STUDY OF AN EXAMPLE: VALIDATION OF A GENERAL CLASS OF PHYTOPLANKTONIC GROWTH MODELS 4.1 Presentation of the model
Theorem 2. (Global qualitative behavior). Almost all trajectories of (~) in a domain (x*) either:
n",.",P
To illustrate the qualitative analysis, we take as an example a class of models describing the behavior of phytoplanktonic cells (whose biomass is X 2) growing on a substrate Xl in a continuous reactor. This bioreactor, for dimensionless variables (for a non zero nutrient supply) can be described by the following system:
n",.",P
• stay in (x*) and go to infinity • stay in n",. ",p(x*) , and goes toward an equilibrium x t in the closure of n",.",p(x*) • go to one of the strict neighbors "'P' (x*) that are accessible.
nu.'
(1)
The variable Xa is the cell quota i.e. the quantity of intracellular nutrient per biomass unit. The functions p and p, represents the absorption rate of the substrate and the growth rate. Among the models ('L.PGM), the Droop model (Droop, 1968; Burmaster, 1979) is largely used in the biological field. For this model we have:
3.1 Graphical representation
We will represent each possible region n",.uP(x*) by a two columns matrix of signs, the first column stands for (jq , and the second for (jp . For example, the region {x E n; Xl > xi , X2 < X 2' Xa > X 3' Xl < 0, X2 < 0, Xa > o} is represented by the matrix:
Hypotheses: To represent growth of phytoplankton, some hypotheses, corroborated by the experiments, are made in the physical domain considered n = {x E R~ ; Xl > 0, X2 > O, xa > a}. • (HI) : The absorption rate p is a non negative bounded function of Xl . It is strictly increasing and verifies p(O) = 0. • (H2) : The growth rate p, is a non negative strictly increasing function of Xs. • (H3): An equilibrium exists in the open domain n.
A possible transition between two regions is represented by an oriented arrow between these regions , the arrows are oriented by the conclusions of theorem (1) . A letter on the arrow will indicate if the variable Xk admits a minimum (mk), a maximum (M k ) , if it crosses its equilibrium xZ top-down (tk) or bottom-up (Tk) ' The set of all n",.up(x*) partitioning connected by the arrows reflecting the transition rules of theorem (1) is called the basic mixed transition graph. The nodes where it is possible to converge to equilibrium will be called equilibrium nodes (simple criteria are given in (Bernard and Gouze, 1997) to identify these nodes).
n
The class of models (~PGM) verifying hypothesis (Hl)-(H2)-(H3) is called the class of Phytoplanktonic Growth Models (PGM). It can easily be verified that the Droop model is in this class. The PGM verify simple properties (see (Bernard and Gouze, 1997) for the proofs) :
Our main theorem has now a "graphical version" (cf. the detailed example below and figure 2) . For the sake of simplicity, we will restrict ourselves to the case where all the trajectories are bounded. We obtain then (cf figure 2):
Property 1. The PGM are loop structured models diagonally x*-monotonous in the domain n (x* is the unique equilibrium in n). Property 2. The trajectories of the PGM are bounded in the positively invariant domain n.
Theorem 3. (Graphical version of the global qualitative theorem, bounded case) On each non equilibrium node, the trajectories follow an arrow of the graph to go toward an accessible node. On each equilibrium node, the trajectories either stay in the node and converge to equilibrium or go toward an accessible node.
4.2 Study of the transition graphs Property 3. There exists 18 feasible qualitative region nu.uP (x*) for the class of PGM . 765
These qualitative situations are the nodes of the graph on fig. l(a), they have been obtained by applying the algorithm of (Bernard and Gouze, 1997). Remark that a priori we have 26 = 64 possible situations, and that we have eliminated 46 cases by this analysis based on the vector field. The consideration of this set of possible domains constitutes therefore a first filter to test the structure of the model. If a qualitative event not belonging to this set is experimentally observed, it means that the system can not rely on the supposed structure. To go further into the description of the qualitative behavior of the model, We can now construct the mixed transition graph by applying Theorem 1 to all the strict neighbors nuqtuPt (x*) and nUq2UP2 (x*) belonging to the set of feasible regions. In the case of this model, let us remark that some strict neighbors have a double simultaneous change of signs: the hyperplane X3 = x3 is included indeed in the nullcline X2 = O.
(a)
1
0 .9 ---..-----1------
Finally we obtain the mixed transition graph (fig. l(a)) associated with the PGM model (I:PGM). This graph summarizes the possible succession of extrema or equilibrium crossings from an initial qualitative situation.
0 .8
---------------------~,:,--...-,.,
.!
8
~
"'., ~
X
~~
o
\d; '''''\. .;). 1, L-"'---_"'---_.a._ .. _..._.' 0-, ' ._ .• ~_:.-:o._---7"0~.:...;'ec<::..;;·.'-"'-"7."-~_'Q<_·'.... ·(~'_·:·.;..,··_o~~. :
·2
·1
0
2
4
5
1lme(days)
(b)
Fig. 1. (a) Basic mixed transition graph. (b) Experiment with the microalgae Dunaliella tertiolecta (with u(O) < 1). The equilibrium X* is experimentally determined (dashed line) .
The PG M dynamics can still be described more precisely by simplifying the transition graph: Property 4. Let u denote the total nutrient concentration in the bioreactor (u = Xl + X2X3). Depending on initial conditions (u(O) < 1 or u(O) > 1), some regions of f! are unreachable, and the transition graph (fig. l(a)) reduces to two simplified graphs (fig. 2).
imum, one maximum, one bottom-up equilibrium crossing and one top-down equilibrium crossing. Qualitative trajectories: Now we can read on the graph the possible succession of qualitative events ; for example, let us start, in the case u(O) < 1 (Fig. 2 up) of the upper left node. The possible sequences are :
Sketch of the proof: Consider the application X - + (x) = u = Xl +X2X3. It is easy to verify that u satisfies a first order differential equation: U 1 - u. The surface Xl + X2X3 1 separates thus the space into 2 unconnected regions + = {x E f!;Xl + X2X3 > 1} and n- = {x E f!;Xl + X2X3 < 1}. Then one has to verify that some domains nUqup (x*) are included in n- or n+. :
=
I
I
0.7
It is noteworthy that the construction of this graph relies only on the sign of the Jacobian matrix, and not on the precise formulation of the model. The biologist can therefore use these graphs as a second filter to test if the basic structure of his model matches the experiments. He can indeed compare the sequences of maxima and minima or the sequences of equilibrium crossings obtained experimentally with the possible sequences contained in the graph.
=
Kz
• (M2, t3), m3, t2, T 1, (m2, T3), M 1, M3 • (M2, t3), t2, m3, T 1, (m2, T3), M I , M3
n
It is thus possible to describe very precisely these trajectories, and to compare them with the data. Let us remark that we have done only qualitative hypotheses on the functions used in the model. The values of the parameters, for example, remain undetermined.
Property 5. The trajectories of the PGM admit at the most for each of the state variables one min-
766
• Finally, some extrema not permitted by the graphs have been observed for biomass . They may also be the consequence of the lag previously described (Cunningham and Maas, 1978). This example illustrates the main advantages of this dynamical qualitative analysis: • It does not need to identify parameters , it is therefore well adapted to biological modeling, where parametric uncertainty is very high . • It enables to compare directly the structure of the model with the data. • It enables therefore to validate not one model, but a class of model, defined by a general structure. • It requires only qualitative features of the data, that can be observed even with noisy data. • It gives necessary conditions on the qualitative behavior of the states variables, that can be useful for validation even when the whole state is not observed.
Acknowledgments : This paper presents research results of the Belgian Programme on Interuniversity Poles of Attraction, initiated by the Belgian State, Prime Minister's Office for Science, Technology and Culture. The scientific responsibility rests with its authors.
The authors acknowledge also the support of the "Environnement, Vie et So cietes' , program and of the GDR 1107 of the C.N.R.S.
Fig. 2. Simplified mixed transition graph for u (O) < 1 (up) u(O) > 1 (down) . 4.3 Qualitative validation
5. REFERENCES Bernard, O. (1995). Etude experimentale et theorique de la croissance de Dunaliella tertiolecta (chlorophyceae) soumise a une limitation variable de nitrate : utilisation de la dynamique transitoire pour la conception et la validation des modeles. PhD thesis. University Pierre et Marie Curie, Paris VI. Bernard, O. and J .-L . Gouze (1995a). Robust validation of uncertain models. In: Proceedings of the Third European Conference on Control (A. Isidori, Ed.). Rome, Italy. pp. 1261-1266. Bernard, O. and J.-L. Gouze (1995b). Transient behavior of biological loop models , with application to the Droop model. Mathematical Biosciences 127, 19-43. Bernard, O. and J.-L. Gouze (1997). Global qualitative behavior of a class of nonlinear biological systems; application to the structural validation of phytoplankton growth models. Rapport de Recherche 3334. INRIA. Bernard, 0 ., G. Malara and A. Sciandra (1996). The effects of a controlled fluctuating nutrient environment on continuous cultures of phyto-
A high number of experiments have been performed with the microalgae Dunaliella terliolecta in a nitrate depleted environment limiting growth (Bernard, 1995) , see figure l(b) . Precise indications on the computer driven environment controlling the culture so as on the conditions of experiments can be found in (Bernard et al., 1996) . From these experiments , three major deficiencies of the Droop model structure have been pointed out. • First, a delay has been observed between extrema of biomass and crossing of the quota equilibrium. We have seen above that , according to the model, these two events should be simultaneous. This lag can be explained by the necessary time for inorganic nitrogen to be incorporated into the organic matter. • Extrema for the substrate have been observed, which are not allowed by the transition graphs. The observed phenomenon may result from the starvation of algae which affects their uptake potential (Collos, 1983) . 767
plankton monitored by a computer. J. Exp. Mar. Biol. Eco1197, 263-278. Burmaster, D. E. (1979). The unsteady continuous culture of phosphate-limited monochrysis Iutheri Droop: experimental and theoretical analysis. J. Exp . Mar. BioI. Ecol. 39(2), 167186. Caswell, H. (1989). Matrix population models. Sinauer Associates Publishers. Collos, Y. (1983) . Transient situations in nitrate assimilation by marine diatoms. 4. Non-linear phenomena and the estimation of the maximum uptake rate. J.Plankton Res. 5, 677-69l. Cunningham, A. and P. Maas (1978) . Time lag and nutrient storage effetcs in the transient growth of Ghlamidomonas reinhardii in nitrogen-limited batch and continuous culture. J. General Microbial. 104, 227-23l. Droop, M. R. (1968) . Vitamin B12 and marine ecology. IV. The kinetics of uptake growth and inhibition in Monochrysis lutheri. J. Mar. Biol. Assoc. 48(3) , 689-733. Goodwin, B. C. (1965) . Oscillatory behaviour in enzymatic control processes. In: Advances in Enzymatic Regulation (G Weber, Ed.). Pergamon. Oxford. Levine, D. S. (1985) . Qualitative theory of a third order nonlinear system with examples in population dynamics and chemical kinetics. Mathematical Biosciences 77, 17-33.
768
QUALITATIVE DYNAMICS OF A CLASS OF NONLINEAR BIOLOGICAL SYSTEMS Olivier BERNARD· Jean-Luc GOUZE ** • GESAME, UGL , Av. G. Lemaitre 4, B-1348 Louvain-La-Neuve, Belgium, [email protected] ** GOMORE Project, INRIA , BP 93, 06902 Sophia-Antipolis, France, [email protected]
Abstract: In this paper we propose a methodology to analyze the global qualitative dynamics of a nonlinear class of differential systems which encompass numerous biological models. We show that, independently of parameters values, the possible succession with respect to time of some qualitative events characterizing transients of state variables (extrema, position ... ) are strongly related to the signs of the Jacobian matrix. The method is illustrated with a set of models usually used to describe phytoplanctonic growth in a bioreactor. The corresponding transition graph is derived and compared with experimental data. This methodology applies well when some parameters are unknown , which is often the case in biological models. Copyright © 1998 IFAC
Resume : Nous proposons une methode d 'analyse de la dynamique qualitative globale d'une classe de systemes differentiels non-lineaires, decrivant de nombreux modeles biologiques. On montre que la succession temporelle de certains evenements dynamiques qualitatifs (passage par un maximum ou un equilibre) depend des signes de la matrice jacobienne. On illustre la methode sur un modele de croissance du phytoplancton dans un bioreacteur. On trace le graphe de transition, qui est compare aux donnees experimentales. Cette methode ne depend pas de la valeur precise des parametres. Keywords: qualitative analysis, nonlinear systems, validation, ecology, biotechnology.
1. INTRODUCTION
systems with monotonous interactions) which has numerous applications in the biological field, including gene regulation models, compartmental systems (Goodwin, 1965j Levine, 1985), cellular growth (Burmaster, 1979), and development of stage structured populations (CasweU, 1989).
In this paper we propose to study the link between the structure of an autonomous differential system (defined by the signs of the Jacobian matrix and its zeros) and its dynamical behavior, in terms of qualitative features for the state variables. The main idea is to describe independently of parameters the dynamics allowed by the structure of a model in order to validate (or invalidate) this structure after a data analysis. In particular, the precise value of some parameters may be unknown, as it is often the case in biological models. We focus more particularly on a class of autonomous differential systems (loop structured
For this class of models, a methodology has recently been proposed (Bemard and Gouze, 1995bj Bemard and Gouze, 1995a) to derive the set of possible qualitative behavior for the state variables. The transient behavior is described either by their possible successive tendencies (increasing or decreasing) (Bemard and Gouze, 1995b) , or by their possible successive position toward
763
O"J E {-I , 1},q E {1 , ... , 2n }}
equilibrium value (larger or lower) (Bemard and Gouze, 1995a). From this analysis it is possible to test experimentally if the extradiagonal signs of the Jacobian matrix are consistent with the data.
Definition 3. For x* E 0 the system with monotonous interactions (E) is diagonally x*-monotonous if for all j in {I , ... , n} , for all q in {I , ... , 211. }, the sign of the partial derivative 8Ii/8xj(x) is fixed on each domain:
In this paper we generalize this result by considering a more precise qualitative state describing simultaneously the tendency of the state and its position toward equilibrium. We consider thus the domain 0 as partitioned by the different regions delimited by the nullclines and the hyperplanes associated to an equilibrium x*. We give in section 2 the theorems that describe the global qualitative behavior of the system from structural considerations i. e. the way trajectories go from one of these regions to another. We present also a graphical version of these theorems, easier to use for validation purposes. Finally, in section 3 we study a class of systems modeling growth of phytoplankton in a continuous bioreactor (chemostat) . The qualitative transient behavior is presented, and some results are derived. It is shown that at the most one maximum, one minimum, one equilibrium crossing bottom-up and one top-down for each state variable is possible. The transition graph obtained is then compared to some experiments, and conclusions on the structure of the models usually used to describe phytoplanktonic growth are drawn.
The domains WO". (x*) are called the orthants of the x*-deviation space or x*-orthants. In the same way, we define the orthants of the velocity space or z-orthants:
Definition 4. For (O"q , O"P) E S~ let us define the following open set: 00"'0"1' (x*) ~f WO" . (x*) n ZO"1'
{x E WO".(x*) ; diag(O"P)I(x)
> O}
Note that some 00"00"1' (x* ) are empty. The remaining non empty 00"' 0"1' (x*) represents therefore the qualitative situations allowed by the model and will be called the possible regions. It is shown in (Bernard and Gouze, 1997) how to derive these possible regions from the signs of the J acobian matrix. We remove from the 00"' 0"1' (x* ) some manifolds of zero measure containing trajectories that may have undesiderable behaviors, and we denote the corresponding non empty domains by nO".O"1'(x*) . We will then show that the transition between these domains obey some rules determined by the extradiagonal terms of the J acobian matrix.
Due to the lack of space, the proofs and details are to be found elsewhere (Bernard and Gouze, 1995b; Bernard and Gouze, 1997) .
2. DEFINITIONS Notations. The notations x > 0 for x = \Xl , ... ,x n ) E Rn means that for all i, Xi> O. For x ERn, sign(x) is the vector with components sign(xi} (with values taken in {-l , O, l}).
Let 0 be an open convex domain of Rn and
=
I a
3. THEOREMS
Cl mapping from 0 onto Rn. We consider on 0 the autonomous differential system :
Theorem 1. (Transitions between regions). Consider a loop system (E) with monotonous interactions and an equilibrium point x*. Suppose that nO"n 0"1'1 (x*) and n".n".1' 2 (x*) are two strict neighbors (i.e. the manifold which separates these regions is of dimension (n-l)) . We denote tk,k+l the sign of the (k, k + 1) element of the Jacobian matrix. A trajectory can go from nO" ' 2".1'2 (x*) to n".n ".1'1 (x*) (resp n""l ".1'1 (x*) to nO"Q2 ".1'2 (x*)] if one of the two following conditions is verified:
(E){x=I(x) Definition 1. A system (E) has a loop structure if hex) = h(Xi , Xi+r) Vi E {1, ... ,n}. Definition 2. The system (E) has monotonous interactions on 0 if each partial derivative 8 Id 8x j (x) for i :f. j never cancels on O.
= =
• tk , k+lO":~l 0":1 (resp. -0":1]; then Xk admits a minimum [resp. a maximum] . • tk,k+l 0"~~1 0"~1 [resp _0"~1]; then Xk crosses its equilibrium xA; bottom-up Iresp top-down].
The signs of the elements define what we call the structure of the system. We will consider the set Sn containing 2n elements:
We say that n""lO"1'l (x*) Iresp accessible from n".'2".1'2 (x*) Iresp 764
n".'2".1'2 nO"n".1'l
(x*)] is (x*)] .
The following theorem describes the behavior of the trajectories of a differential system (~), loop structured and diagonally-monotonous, with monotonous interactions.
4. STUDY OF AN EXAMPLE: VALIDATION OF A GENERAL CLASS OF PHYTOPLANKTONIC GROWTH MODELS 4.1 Presentation of the model
Theorem 2. (Global qualitative behavior). Almost all trajectories of (~) in a domain (x*) either:
n",.",P
To illustrate the qualitative analysis, we take as an example a class of models describing the behavior of phytoplanktonic cells (whose biomass is X 2) growing on a substrate Xl in a continuous reactor. This bioreactor, for dimensionless variables (for a non zero nutrient supply) can be described by the following system:
n",.",P
• stay in (x*) and go to infinity • stay in n",. ",p(x*) , and goes toward an equilibrium x t in the closure of n",.",p(x*) • go to one of the strict neighbors "'P' (x*) that are accessible.
nu.'
(1)
The variable Xa is the cell quota i.e. the quantity of intracellular nutrient per biomass unit. The functions p and p, represents the absorption rate of the substrate and the growth rate. Among the models ('L.PGM), the Droop model (Droop, 1968; Burmaster, 1979) is largely used in the biological field. For this model we have:
3.1 Graphical representation
We will represent each possible region n",.uP(x*) by a two columns matrix of signs, the first column stands for (jq , and the second for (jp . For example, the region {x E n; Xl > xi , X2 < X 2' Xa > X 3' Xl < 0, X2 < 0, Xa > o} is represented by the matrix:
Hypotheses: To represent growth of phytoplankton, some hypotheses, corroborated by the experiments, are made in the physical domain considered n = {x E R~ ; Xl > 0, X2 > O, xa > a}. • (HI) : The absorption rate p is a non negative bounded function of Xl . It is strictly increasing and verifies p(O) = 0. • (H2) : The growth rate p, is a non negative strictly increasing function of Xs. • (H3): An equilibrium exists in the open domain n.
A possible transition between two regions is represented by an oriented arrow between these regions , the arrows are oriented by the conclusions of theorem (1) . A letter on the arrow will indicate if the variable Xk admits a minimum (mk), a maximum (M k ) , if it crosses its equilibrium xZ top-down (tk) or bottom-up (Tk) ' The set of all n",.up(x*) partitioning connected by the arrows reflecting the transition rules of theorem (1) is called the basic mixed transition graph. The nodes where it is possible to converge to equilibrium will be called equilibrium nodes (simple criteria are given in (Bernard and Gouze, 1997) to identify these nodes).
n
The class of models (~PGM) verifying hypothesis (Hl)-(H2)-(H3) is called the class of Phytoplanktonic Growth Models (PGM). It can easily be verified that the Droop model is in this class. The PGM verify simple properties (see (Bernard and Gouze, 1997) for the proofs) :
Our main theorem has now a "graphical version" (cf. the detailed example below and figure 2) . For the sake of simplicity, we will restrict ourselves to the case where all the trajectories are bounded. We obtain then (cf figure 2):
Property 1. The PGM are loop structured models diagonally x*-monotonous in the domain n (x* is the unique equilibrium in n). Property 2. The trajectories of the PGM are bounded in the positively invariant domain n.
Theorem 3. (Graphical version of the global qualitative theorem, bounded case) On each non equilibrium node, the trajectories follow an arrow of the graph to go toward an accessible node. On each equilibrium node, the trajectories either stay in the node and converge to equilibrium or go toward an accessible node.
4.2 Study of the transition graphs Property 3. There exists 18 feasible qualitative region nu.uP (x*) for the class of PGM . 765
These qualitative situations are the nodes of the graph on fig. l(a), they have been obtained by applying the algorithm of (Bernard and Gouze, 1997). Remark that a priori we have 26 = 64 possible situations, and that we have eliminated 46 cases by this analysis based on the vector field. The consideration of this set of possible domains constitutes therefore a first filter to test the structure of the model. If a qualitative event not belonging to this set is experimentally observed, it means that the system can not rely on the supposed structure. To go further into the description of the qualitative behavior of the model, We can now construct the mixed transition graph by applying Theorem 1 to all the strict neighbors nuqtuPt (x*) and nUq2UP2 (x*) belonging to the set of feasible regions. In the case of this model, let us remark that some strict neighbors have a double simultaneous change of signs: the hyperplane X3 = x3 is included indeed in the nullcline X2 = O.
(a)
1
0 .9 ---..-----1------
Finally we obtain the mixed transition graph (fig. l(a)) associated with the PGM model (I:PGM). This graph summarizes the possible succession of extrema or equilibrium crossings from an initial qualitative situation.
0 .8
---------------------~,:,--...-,.,
.!
8
~
"'., ~
X
~~
o
\d; '''''\. .;). 1, L-"'---_"'---_.a._ .. _..._.' 0-, ' ._ .• ~_:.-:o._---7"0~.:...;'ec<::..;;·.'-"'-"7."-~_'Q<_·'.... ·(~'_·:·.;..,··_o~~. :
·2
·1
0
2
4
5
1lme(days)
(b)
Fig. 1. (a) Basic mixed transition graph. (b) Experiment with the microalgae Dunaliella tertiolecta (with u(O) < 1). The equilibrium X* is experimentally determined (dashed line) .
The PG M dynamics can still be described more precisely by simplifying the transition graph: Property 4. Let u denote the total nutrient concentration in the bioreactor (u = Xl + X2X3). Depending on initial conditions (u(O) < 1 or u(O) > 1), some regions of f! are unreachable, and the transition graph (fig. l(a)) reduces to two simplified graphs (fig. 2).
imum, one maximum, one bottom-up equilibrium crossing and one top-down equilibrium crossing. Qualitative trajectories: Now we can read on the graph the possible succession of qualitative events ; for example, let us start, in the case u(O) < 1 (Fig. 2 up) of the upper left node. The possible sequences are :
Sketch of the proof: Consider the application X - + (x) = u = Xl +X2X3. It is easy to verify that u satisfies a first order differential equation: U 1 - u. The surface Xl + X2X3 1 separates thus the space into 2 unconnected regions + = {x E f!;Xl + X2X3 > 1} and n- = {x E f!;Xl + X2X3 < 1}. Then one has to verify that some domains nUqup (x*) are included in n- or n+. :
=
I
I
0.7
It is noteworthy that the construction of this graph relies only on the sign of the Jacobian matrix, and not on the precise formulation of the model. The biologist can therefore use these graphs as a second filter to test if the basic structure of his model matches the experiments. He can indeed compare the sequences of maxima and minima or the sequences of equilibrium crossings obtained experimentally with the possible sequences contained in the graph.
=
Kz
• (M2, t3), m3, t2, T 1, (m2, T3), M 1, M3 • (M2, t3), t2, m3, T 1, (m2, T3), M I , M3
n
It is thus possible to describe very precisely these trajectories, and to compare them with the data. Let us remark that we have done only qualitative hypotheses on the functions used in the model. The values of the parameters, for example, remain undetermined.
Property 5. The trajectories of the PGM admit at the most for each of the state variables one min-
766
• Finally, some extrema not permitted by the graphs have been observed for biomass . They may also be the consequence of the lag previously described (Cunningham and Maas, 1978). This example illustrates the main advantages of this dynamical qualitative analysis: • It does not need to identify parameters , it is therefore well adapted to biological modeling, where parametric uncertainty is very high . • It enables to compare directly the structure of the model with the data. • It enables therefore to validate not one model, but a class of model, defined by a general structure. • It requires only qualitative features of the data, that can be observed even with noisy data. • It gives necessary conditions on the qualitative behavior of the states variables, that can be useful for validation even when the whole state is not observed.
Acknowledgments : This paper presents research results of the Belgian Programme on Interuniversity Poles of Attraction, initiated by the Belgian State, Prime Minister's Office for Science, Technology and Culture. The scientific responsibility rests with its authors.
The authors acknowledge also the support of the "Environnement, Vie et So cietes' , program and of the GDR 1107 of the C.N.R.S.
Fig. 2. Simplified mixed transition graph for u (O) < 1 (up) u(O) > 1 (down) . 4.3 Qualitative validation
5. REFERENCES Bernard, O. (1995). Etude experimentale et theorique de la croissance de Dunaliella tertiolecta (chlorophyceae) soumise a une limitation variable de nitrate : utilisation de la dynamique transitoire pour la conception et la validation des modeles. PhD thesis. University Pierre et Marie Curie, Paris VI. Bernard, O. and J .-L . Gouze (1995a). Robust validation of uncertain models. In: Proceedings of the Third European Conference on Control (A. Isidori, Ed.). Rome, Italy. pp. 1261-1266. Bernard, O. and J.-L. Gouze (1995b). Transient behavior of biological loop models , with application to the Droop model. Mathematical Biosciences 127, 19-43. Bernard, O. and J.-L. Gouze (1997). Global qualitative behavior of a class of nonlinear biological systems; application to the structural validation of phytoplankton growth models. Rapport de Recherche 3334. INRIA. Bernard, 0 ., G. Malara and A. Sciandra (1996). The effects of a controlled fluctuating nutrient environment on continuous cultures of phyto-
A high number of experiments have been performed with the microalgae Dunaliella terliolecta in a nitrate depleted environment limiting growth (Bernard, 1995) , see figure l(b) . Precise indications on the computer driven environment controlling the culture so as on the conditions of experiments can be found in (Bernard et al., 1996) . From these experiments , three major deficiencies of the Droop model structure have been pointed out. • First, a delay has been observed between extrema of biomass and crossing of the quota equilibrium. We have seen above that , according to the model, these two events should be simultaneous. This lag can be explained by the necessary time for inorganic nitrogen to be incorporated into the organic matter. • Extrema for the substrate have been observed, which are not allowed by the transition graphs. The observed phenomenon may result from the starvation of algae which affects their uptake potential (Collos, 1983) . 767
plankton monitored by a computer. J. Exp. Mar. Biol. Eco1197, 263-278. Burmaster, D. E. (1979). The unsteady continuous culture of phosphate-limited monochrysis Iutheri Droop: experimental and theoretical analysis. J. Exp . Mar. BioI. Ecol. 39(2), 167186. Caswell, H. (1989). Matrix population models. Sinauer Associates Publishers. Collos, Y. (1983) . Transient situations in nitrate assimilation by marine diatoms. 4. Non-linear phenomena and the estimation of the maximum uptake rate. J.Plankton Res. 5, 677-69l. Cunningham, A. and P. Maas (1978) . Time lag and nutrient storage effetcs in the transient growth of Ghlamidomonas reinhardii in nitrogen-limited batch and continuous culture. J. General Microbial. 104, 227-23l. Droop, M. R. (1968) . Vitamin B12 and marine ecology. IV. The kinetics of uptake growth and inhibition in Monochrysis lutheri. J. Mar. Biol. Assoc. 48(3) , 689-733. Goodwin, B. C. (1965) . Oscillatory behaviour in enzymatic control processes. In: Advances in Enzymatic Regulation (G Weber, Ed.). Pergamon. Oxford. Levine, D. S. (1985) . Qualitative theory of a third order nonlinear system with examples in population dynamics and chemical kinetics. Mathematical Biosciences 77, 17-33.
768