3 Modélisation du chémostat

Nous allons maintenant nous intéresser au fonctionnement du chémostat (Fig. 4), en particulier à la modélisation simultanée de la croissance bactérienne (variable \(N(t)\)) et de la consommation en substrat (variable \(S(t)\)) à l’intérieur de la chambre de culture.

Fig. 4 : Schéma de fonctionnement d’un chémostat. Une substance nutritive pénètre dans la chambre de culture bactérienne à la concentration \(S_0\) avec un flux d’entrée \(F\), égal au taux de sortie, de sorte que le volume total \(V\) reste constant. D’après Edelstein-Keshet (1984).

Dans un chémostat, l’objectif est double :

1. Choisir un flux d’entrée de sorte à ne pas lessiver complètement la culture bactérienne et du coup éliminer les bactéries du système ; typiquement, on cherchera à maintenir le système dans un état d’équilibre.

2. Faire en sorte que le réapprovisionnement en nutriments soit suffisamment rapide pour que la croissance bactérienne se fasse normalement.

Ainsi, il faut pouvoir fixer à bon escient les paramètres \(F\), \(S_0\) et \(V\).

3.1 Variables et paramètres

Les variables et les paramètres qui vont nous intéresser dans cette partie sont résumés dans le tableau ci-dessous (Tab. 1).

Tab. 1 : Variables et paramètres du modèle de chémostat.

Quantité	Symbole	Dimension
Concentration en nutriments à l’instant \(t\) dans la chambre de culture	\(S(t)\)	masse/volume
Biomasse bactérienne à l’instant \(t\) dans la chambre de culture	\(N(t)\)	CFU/volume
Concentration en nutriments dans le réservoir d’entrée	\(S_0\)	masse/volume
Volume de la chambre de culture	\(V\)	volume
Flux d’entrée/sortie	\(F\)	volume/temps
Taux de croissance bactérien	\(r\)	1/temps
Constante de rendement (cf. § 3.2)	\(Y=1/\alpha\)	masse/volume

3.2 Hypothèses et premières équations

Un certain nombre d’hypothèses sont par ailleurs nécessaires à la construction du modèle :

1. Le milieu de culture est infiniment mélangé et donc homogène. Ceci nous autorise à utiliser pour la modélisation des équations différentielles ordinaires avec le temps comme seule variable indépendante.

On peut dès lors proposer une première équation pour décrire l’évolution de la densité bactérienne :

\[\frac{dN(t)}{dt} = r N(t) - F \frac{N(t)}{V}\] La partie \(r N(t)\) correspond à la croissance exponentielle par unité de temps avec \(r\) (dont l’unité est 1/temps) le taux de croissance malthusien (§ 2.2).

La partie \(- F \frac{N(t)}{V}\) correspond à ce qui sort de la cuve par unité de temps.

La division par \(V\) assure le respect des dimensions dans l’équation. À cause de l’unité de \(F\), la quantité \(\frac{F}{V}\) s’exprime en \([temps]^{-1}\).

2. Par soucis de simplification, on considère en première approximation qu’il n’y a qu’un seul type de nutriments dans la chambre de culture, dont la concentration influence la croissance bactérienne.

3. Le taux de croissance bactérien dépend de la disponibilité en nutriments (Monod, 1942) de sorte que :

\[r = r(S)\]

Nous verrons ultérieurement quelle hypothèse choisir pour la fonction \(r(S)\) : une simple proportionnalité ou un modèle plus complexe.

Il nous faut maintenant une deuxième équation pour décrire l’évolution au cours du temps de la concentration en nutriments.

4. Nous allons supposer que la consommation en nutriments se fait en continu du fait de la croissance bactérienne, et que \(\alpha\) unités de nutriment sont consommées pour produire 1 unité de bactéries. La quantité \(Y=1/\alpha\) représente ainsi le rendement, et on peut écrire :

\[\frac{dS(t)}{dt} = -\alpha r(S) N(t) - F \frac{S(t)}{V} + F \frac{S_0}{V}\]

avec \(-\alpha r(S) N(t)\) ce qui est consommé, \(- F \frac{S(t)}{V}\) ce qui sort de la cuve et \(+ F \frac{S_0}{V}\) ce qui entre dans la cuve, par unité de temps.

3.3 La notion de substrat limitant

En milieu confiné, la croissance bactérienne peut se trouver retardée et/ou limitée (Monod, 1935) et ce pour différentes raisons : le substrat peut s’épuiser, il peut y avoir accumulation de toxines produites par les bactéries elles-mêmes, des variations de pH peuvent être induites par un produit de fermentation (e.g., acide lactique, acide acétique), ou bien encore l’oxygène peut s’épuiser (Lobry, 1991).

Dans les conditions qui prévalent en milieu liquide synthétique (comme dans un chémostat), c’est l’épuisement du substrat limitant qui arrête la croissance bactérienne. Ainsi, Monod (Monod, 1935 ; Monod et Teissier, 1936) proposa de supposer que la variation du taux de croissance bactérien (ici \(r(S)\)) en fonction de la concentration en nutriments (\(S\)) n’est pas proportionnel à \(S\) mais suit une relation monotone croissante :

\[r(S) = \frac{r_{max} S}{K_s + S}\]

Cette équation ressemble à s’y méprendre à celle du modèle de Michaelis-Menten (ou de Michaelis-Menten-Henri) permettant de décrire la cinétique d’une réaction chimique catalysée par une enzyme agissant sur un substrat unique pour donner un produit. Le modèle de Michaelis-Menten relie la vitesse de la réaction chimique \(v\) à la concentration en substrat \([S]\) ; il inclut des paramètres constants caractéristiques de l’enzyme :

\[v([S]) = \frac{v_{max} [S]}{K_M + [S]}\]

\(K_M\) est la constante de Michaelis spécifique de l’enzyme ; elle a la dimension d’une concentration (même unité que \([S]\)). Elle est définie par \(K_M = \frac{k_3 + k_2}{k_1}\) où les constantes \(k_i\) sont les constantes de vitesse de la réaction chimique suivante :

\[E+S \, \stackrel {k_1}{\underset {k_2}{\rightleftharpoons}} \, ES \, \stackrel{k_3}{\rightarrow} \, E + P\] \[v = k_3 \times [ES]\]

\(K_M\) est aussi la concentration en substrat pour laquelle la vitesse de la réaction chimique est égale à la moitié de la vitesse maximale.

• Ecrire le code R permettant de reproduire la figure ci-dessous (Fig. 5), que vous compléterez.

Fig. 5 : Représentation graphique du taux de croissance bactérien (\(r\)) en fonction de la concentration en nutriments (\(S\)).

• A quoi correspondent mathématiquement les paramètres \(r_{max}\) et \(K_s\) ?

• Quelle est la signification biologique du paramètre \(r_{max}\) ?

3.4 Reparamétrisation du modèle

Sur la base des deux équations précédentes, on obtient finalement un système de deux EDO couplées pour décrire le fonctionnement du chémostat : \[\left\{ \begin{array}{l} \frac{{dN\left( t \right)}}{{dt}} = \left( {\frac{{{r_{\max }}S\left( t \right)}}{{{K_s} + S\left( t \right)}}} \right)N\left( t \right) - \frac{F}{V}N\left( t \right)\\ \frac{{dS\left( t \right)}}{{dt}} = - \alpha \left( {\frac{{{r_{\max }}S\left( t \right)}}{{{K_s} + S\left( t \right)}}} \right)N\left( t \right) - \frac{F}{V}S\left( t \right) + \frac{F}{V}{S_0} \end{array} \right.\] Ce système contient 6 paramètres : \(r_{max}\), \(K_s\), \(F\), \(V\), \(\alpha\) et \(S_0\). Pour simplifier son analyse, nous allons tenter de trouver une nouvelle paramétrisation par un changement des variables sans dimension.

Supposons que l’on puisse écrire \(N(t) = N^*(t) \times \hat N\), \(S(t) = S^*(t) \times \hat S\) et \(t = t^* \times \tau\), avec \(N^*(t)\), \(S^*(t)\) et \(t^*\) de nouvelles variables sans dimension, et \(\hat N\), \(\hat S\) et \(\tau\) des constantes indépendantes du temps. Ainsi : \[\left\{ \begin{array}{l} \frac{{d\left( {{N^*}\hat N} \right)}}{{d\left( {{t^*}\tau } \right)}} = \left( {\frac{{{r_{\max }}{S^*}\hat S}}{{{K_s} + {S^*}\hat S}}} \right){N^*}\hat N - \frac{F}{V}{N^*}\hat N\\ \frac{{d\left( {{S^*}\hat S} \right)}}{{d\left( {{t^*}\tau } \right)}} = - \alpha \left( {\frac{{{r_{\max }}{S^*}\hat S}}{{{K_s} + {S^*}\hat S}}} \right){N^*}\hat N - \frac{F}{V}{S^*}\hat S + \frac{F}{V}{S_0} \end{array} \right.\]

Pour comprendre ce que sont ces nouvelles variables, prenons l’exemple d’une biomasse bactérienne qui vaut \(N = 10^5\) UFC/ml, alors : \[\begin{array}{l} N = {10^5} \quad {\rm{ UFC/ml}}\\ = 1\;\left( { \times \;{{10}^5}\;{\rm{UFC}}} \right){\rm{/ml}}\\ = {\rm{100}}\quad \,{\rm{UFC/}}\mu {\rm{l}}\\ = {N^*} \times \hat N \end{array}\] avec \(N^*\) un nombre sans dimension et \(\hat N\) une quantité qui caractérise l’unité de mesure selon des échelles différentes.

Dans le système précédent, si, de part et d’autre du signe \(=\), on multiplie par \(\tau\) et on divise par \(\hat N\) dans la première équation (resp. par \(\hat S\) dans la deuxième équation), on obtient : \[\left\{ \begin{array}{l} \frac{{d\left( {{N^*}} \right)}}{{d\left( {{t^*}} \right)}} = \tau {r_{\max }}\left( {\frac{{{S^*}}}{{{K_s}/\hat S + {S^*}}}} \right){N^*} - \tau \frac{F}{V}{N^*}\\ \frac{{d\left( {{S^*}} \right)}}{{d\left( {{t^*}} \right)}} = - \alpha \tau {r_{\max }}\frac{{\hat N}}{{\hat S}}\left( {\frac{{{S^*}}}{{{K_s}/\hat S + {S^*}}}} \right){N^*} - \tau \frac{F}{V}{S^*} + \tau \frac{F}{V}\frac{{{S_0}}}{{\hat S}} \end{array} \right.\]

Si on pose maintenant \(\tau = \frac{V}{F}\), \(\hat S = K_s\) et \(\hat N = \frac{K_s}{\alpha \tau r_{max}}\), alors il vient : \[\left\{ \begin{array}{l} \frac{{d{N^*}}}{{d{t^*}}} = \tau {r_{\max }}\left( {\frac{{{S^*}}}{{1 + {S^*}}}} \right){N^*} - {N^*}\\ \frac{{d{S^*}}}{{d{t^*}}} = - \left( {\frac{{{S^*}}}{{1 + {S^*}}}} \right){N^*} - {S^*} + \frac{{{S_0}}}{{K_s}} \end{array} \right.\]

• Quelle interprétation peut-on donner au paramètre \(\tau = \frac{V}{F}\) ?

Enfin, si on pose \(\color{green}{\alpha_1 = \tau r_{max} = \frac{V r_{max}}{F}}\) et \(\color{green}{\alpha_2 = \frac{S_0}{K_s}}\), alors on obtient un nouveau modèle qui ne contient plus que deux paramètres : \[\left\{ \begin{array}{l} \frac{{d{N^*}}}{{d{t^*}}} = {\color{green}{\alpha _1}}\left( {\frac{{{S^*}}}{{1 + {S^*}}}} \right){N^*} - {N^*}\\ \frac{{d{S^*}}}{{d{t^*}}} = - \left( {\frac{{{S^*}}}{{1 + {S^*}}}} \right){N^*} - {S^*} + {\color{green}{\alpha _2}} \end{array} \right.\]

C’est ce modèle que nous allons maintenant étudier, d’abord par une analyse mathématique, puis ensuite avec des simulations numériques. Par commodité pour la suite, nous enlèverons les \(^*\) pour désigner les variables.

• Quelle est la signification du paramètre \(\alpha_2\) ?

3.5 Analyse mathématique du modèle du chémostat

3.5.1 Recherche des points d’équilibre

Conformément à ce qui a été dit plus haut, l’objectif est d’arriver à maintenir le chémostat en état d’équilibre. Cet état d’équilibre est atteint si les variables du système, c’est-à-dire \(N(t)\) et \(S(t)\) restent constantes toutes les deux simultanément. Mathématiquement, cela revient donc à vérifier :

\[\left\{ \begin{array}{l} \frac{{dN}}{{dt}} = 0\\ \frac{{dS}}{{dt}} = 0 \end{array} \right. \Leftrightarrow \left\{ \begin{array}{l} {\alpha _1}\left( {\frac{S}{{1 + S}}} \right)N - N = 0\\ - \left( {\frac{S}{{1 + S}}} \right)N - S + {\alpha _2} = 0 \end{array} \right.\]

On qualifie de points d’équilibre les solutions du système \(\frac{dN}{dt} = 0\) et \(\frac{dS}{dt} = 0\).

La première équation conduit aux deux solutions \(N_1=0\) ou \(\frac{S_2}{1+S_2}=\frac{1}{\alpha_1}\), c’est-à-dire \(N_1=0\) et \(S_2 = \frac{1}{\alpha_1 - 1}\).

On considérera dans la suite que \(\alpha_1 \neq 1\) pour que \(S_2\) existe mathématiquement.

• Si \(N_1=0\), alors la deuxième équation du système conduit à \(S_1 = \alpha_2\). Dans le plan \((N, S)\), on a donc un premier point d’équilibre \(\color{green}{(N_1, S_1) = (0, \alpha_2)}\).

  • Quelle interprétation biologique pouvez-vous donner à ce point d’équilibre ?

• Si \(S_2 = \frac{1}{\alpha_1 - 1}\), alors la deuxième équation du système conduit à \(\left(\frac{S_2}{1+S_2}\right)N_2 = \alpha_2 - S_2\), c’est-à-dire \(N_2 = \left(\frac{1+S_2}{S_2}\right)(\alpha_2 - S_2) = \alpha_1(\alpha_2 - S_2)\). On obtient donc un deuxième point d’équilibre \(\color{green}{(N_2, S_2) = (\alpha_1(\alpha_2-\frac{1}{\alpha_1 - 1}), \frac{1}{\alpha_1 - 1})}\).

  • A quelle(s) condition(s) sur les paramètres ce deuxième point d’équilibre a-t-il du sens biologiquement ?

  • Quelle interprétation biologique donnez-vous à ce deuxième point d’équilibre ?

3.5.2 Stabilité des points d’équilibre

Maintenant que les points d’équilibre ont été identifiés, il est intéressant de s’intéresser à leur stabilité. Si pour une raison ou une autre, l’état d’équilibre est perturbé, le système y revient-il rapidement ou bien la dynamique du système change-t-elle radicalement ?

La théorie des systèmes dynamiques stipule qu’un point d’équilibre \((x_{eq}, y_{eq})\) du système d’équations : \[\left\{ \begin{array}{l} \frac{dx}{dt} = f(x,y)\\ \frac{dy}{dt} = g(x,y)\\ \end{array} \right.\] est asymptotiquement stable si et seulement si les deux conditions suivantes sont réunies :

\[\left\{ {\begin{array}{*{20}{l}} {\frac{{\partial f({x_{eq}},{y_{eq}})}}{{\partial x}} + \frac{{\partial g({x_{eq}},{y_{eq}})}}{{\partial y}} < 0 \quad \color{magenta}{({\rm{critère-1)}}}}\\ {\frac{{\partial f({x_{eq}},{y_{eq}})}}{{\partial x}}.\frac{{\partial g({x_{eq}},{y_{eq}})}}{{\partial y}} - \frac{{\partial f({x_{eq}},{y_{eq}})}}{{\partial y}}.\frac{{\partial g({x_{eq}},{y_{eq}})}}{{\partial x}} > 0 \quad \color{magenta}{({\rm{critère-2)}}}} \end{array}} \right.\]

où les termes \(\frac{\partial f(x_{eq},y_{eq})}{\partial x}\), \(\frac{\partial f(x_{eq},y_{eq})}{\partial y}\), \(\frac{\partial g(x_{eq},y_{eq})}{\partial x}\), et \(\frac{\partial g(x_{eq},y_{eq})}{\partial y}\) sont les dérivées partielles des fonctions \(f(x,y)\) et \(g(x,y)\) par rapport aux variables \(x\) et \(y\), respectivement, évaluées en \(x=x_{eq}\) et \(y=y_{eq}\).

Dans le cas du modèle du chémostat, nous avons : \[\begin{array}{l} f\left( {N,S} \right) = {\alpha _1}\left( {\frac{S}{{1 + S}}} \right)N - N\\ g\left( {N,S} \right) = - \left( {\frac{S}{{1 + S}}} \right)N - S + {\alpha _2} \end{array}\]

• Vérifiez que \(\left(\frac{x}{1+x} \right)' = \frac{1}{(1+x)^2}\)

Pour le modèle du chémostat, le calcul des dérivées partielles conduit à : \[\begin{array}{l} \frac{{\partial f\left( {N,S} \right)}}{{\partial N}} = {\alpha _1}\left( {\frac{S}{{1 + S}}} \right) - 1\\ \frac{{\partial f\left( {N,S} \right)}}{{\partial S}} = {\alpha _1}\frac{N}{{{{\left( {1 + S} \right)}^2}}}\\ \frac{{\partial g\left( {N,S} \right)}}{{\partial N}} = - \left( {\frac{S}{{1 + S}}} \right)\\ \frac{{\partial g\left( {N,S} \right)}}{{\partial S}} = - \frac{N}{{{{\left( {1 + S} \right)}^2}}} - 1 \end{array}\]

Il nous faut maintenant calculer ces dérivés partielles en chacun des points d’équilibre \((N_1,S_1)\) et \((N_2,S_2)\). Pour simplifier les calculs, posons : \[A = \frac{\alpha_2}{1+\alpha_2}\] et \[B = \frac{N_2}{(1+S_2)^2}\] Les deux quantités \(A\) et \(B\) sont \(>0\).

• Au point d’équilibre \((N_1, S_1) = (0, \alpha_2)\), il vient : \[\begin{array}{l} \frac{{\partial f\left( {{N_1},{S_1}} \right)}}{{\partial N}} = {\alpha _1}\left( {\frac{{{\alpha _2}}}{{1 + {\alpha _2}}}} \right) - 1 = {\alpha _1}A - 1\\ \frac{{\partial f\left( {{N_1},{S_1}} \right)}}{{\partial S}} = 0\\ \frac{{\partial g\left( {{N_1},{S_1}} \right)}}{{\partial N}} = - \left( {\frac{{{\alpha _2}}}{{1 + {\alpha _2}}}} \right) = - A\\ \frac{{\partial g\left( {{N_1},{S_1}} \right)}}{{\partial S}} = - 1 \end{array}\]

Ainsi : \[\begin{array}{l} \frac{{\partial f\left( {{N_1},{S_1}} \right)}}{{\partial N}} + \frac{{\partial g\left( {{N_1},{S_1}} \right)}}{{\partial S}} = {\alpha _1}A - 2\\ \frac{{\partial f\left( {{N_1},{S_1}} \right)}}{{\partial N}}.\frac{{\partial g\left( {{N_1},{S_1}} \right)}}{{\partial S}} - \frac{{\partial f\left( {{N_1},{S_1}} \right)}}{{\partial S}}.\frac{{\partial g\left( {{N_1},{S_1}} \right)}}{{\partial N}} = 1 - {\alpha _1}A \end{array}\]

Comme nous l’avons vu précédemment, le deuxième point d’équilibre a du sens biologiquement si \(\alpha_1 > 1\) et \(\alpha_2(\alpha_1 - 1) > 1\). Nous voyons donc que si ces deux conditions sont vérifiées, alors : \[1 - {\alpha _1}A = 1 - {\alpha _1}\left( {\frac{{{\alpha _2}}}{{1 + {\alpha _2}}}} \right) = \frac{{1 + {\alpha _2} - {\alpha _1}{\alpha _2}}}{{1 + {\alpha _2}}} = \frac{{1 - {\alpha _2}\left( {{\alpha _1} - 1} \right)}}{{1 + {\alpha _2}}} < 0\] Par conséquent, le point d’équilibre \((N_1, S_1)\) est instable.

• Le deuxième point d’équilibre \((N_2, S_2)\) s’écrit \((\alpha_1(\alpha_2-\frac{1}{\alpha_1 - 1}), \frac{1}{\alpha_1 - 1})\). Il n’a de sens biologique que si \(\alpha_1 > 1\) et \(\alpha_2(\alpha_1 - 1) > 1\), et il vérifie \(\frac{S_2}{1+S_2} = \frac{1}{\alpha_1}\). Ainsi, il vient : \[\begin{array}{l} \frac{{\partial f\left( {{N_2},{S_2}} \right)}}{{\partial N}} = 0\\ \frac{{\partial f\left( {{N_2},{S_2}} \right)}}{{\partial S}} = {\alpha _1}B\\ \frac{{\partial g\left( {{N_2},{S_2}} \right)}}{{\partial N}} = - \frac{1}{{{\alpha _1}}}\\ \frac{{\partial g\left( {{N_2},{S_2}} \right)}}{{\partial S}} = - B - 1 \end{array}\]

Le calcul des conditions de stabilité conduit à : \[\begin{array}{l} \frac{{\partial f\left( {{N_2},{S_2}} \right)}}{{\partial N}} + \frac{{\partial g\left( {{N_2},{S_2}} \right)}}{{\partial S}} = - B - 1 < 0\\ \frac{{\partial f\left( {{N_2},{S_2}} \right)}}{{\partial N}}.\frac{{\partial g\left( {{N_2},{S_2}} \right)}}{{\partial S}} - \frac{{\partial f\left( {{N_2},{S_2}} \right)}}{{\partial S}}.\frac{{\partial g\left( {{N_2},{S_2}} \right)}}{{\partial N}} = B > 0 \end{array}\] Par conséquent, le point d’équilibre \((N_2, S_2)\) est asymptotiquement stable.

• Que pouvez-vous conclure de cette analyse de stabilité sur la dynamique du chémostat telle que décrite par ce modèle ?

3.5.3 Interprétation des conditions sur les paramètres

Nous avons établi précédemment que le modèle du chémostat présente deux points d’équilibre dont le point d’équilibre \((N_2, S_2)\) qui fait sens biologiquement si \(\alpha_1 > 1\) et \(\alpha_2(\alpha_1 - 1) > 1\) et qui est asymptotiquement stable.

• La condition \(\alpha_1 > 1\) est équivalente à \(r_{max} > \frac{1}{\tau} \Leftrightarrow \frac{1}{r_{max}} < \frac{V}{F}\). Nous avons établi précédemment, dans le cadre du modèle de croissance exponentielle, que le temps de doublement de la population bactérienne était \(\tau_d = \frac{\ln 2}{r_{max}}\). Ainsi, la condition \(\alpha_1 > 1\) équivaut à \(\tau_d < \ln 2 \frac{ V}{F}\), c’est-à-dire qu’elle impose que \(\tau_d\) soit inférieur au temps nécessaire pour vider la chambre de culture (\(\times \ln 2\)), autrement dit les bactéries doivent pouvoir au moins doubler leur population pendant une période de temps égale au temps de vidange de la chambre.

• La seconde condition \(\alpha_2(\alpha_1 - 1) > 1\) équivaut à \(\frac{S_0}{K_s} > \frac{1}{\alpha_1 - 1}\), c’est-à-dire \(\frac{S_0}{K_s} > S_2^*\), soit encore \(S_0 > K_s S_2^*\). Si on se rappelle que \(\hat S = K_s\), il vient \(S_0 > \hat S S_2^*\). La quantité \(\hat S S_2^*\) correspond à la valeur de \(S_2\) exprimée dans l’unité de départ avant le re-dimensionnement, c’est-à-dire exprimée en masse/volume, comme l’est \(S_0\). En conclusion, nous tombons sur un résultat trivial : la concentration en nutriments à l’équilibre (\(\hat S S_2^*\)) ne peut pas être supérieure à la concentration du stock de nutriment (\(S_0\)).

3.6 Représentations graphiques des solutions d’un système dynamique

3.6.1 Le portrait de phase

Dans un système tel que celui que nous avons \[\left\{ \begin{array}{l} \frac{dN}{dt} = f(N, S)\\ \frac{dS}{dt} = g(N, S)\\ \end{array} \right.\] les solutions sont de la forme \((N(t), S(t))\).

Ces solutions peuvent être représentées graphiquement soit dans un plan \((N, S)\), qu’on appelle le plan de phase, soit dans des plans \((t, N(t))\) ou \((t, S(t))\). Dans le plan de phase, la représentation graphique d’une solution \((N(t), S(t))\) dessine une trajectoire depuis une condition initiale \((N(0), S(0))\) et par exemple jusqu’à un point d’équilibre asymptotiquement stable \((N_{eq}, S_{eq})\) (Fig. 6).

Fig. 6 : Représentation graphique d’une solution \((N(t),S(t))\) dans le plan de phase depuis une condition initiale \((N(0),S(0))\) jusqu’à un point d’équilibre asymptotiquement stable \((N_{eq},S_{eq})\).

Sans résoudre analytiquement le système dynamique, la solution \((N(t), S(t))\) est en fait connue indirectement par les fonctions \(f(N, S)\) et \(g(N, S)\) qui caractérisent en chaque point (et donc à chaque pas de temps) le vecteur vitesse : \[\textrm{V} = \left( {\begin{array}{*{20}{c}} {\frac{dN}{dt}}\\ {\frac{dS}{dt}} \end{array}} \right)\] qui est tangent à la trajectoire (cf. Fig. 6, point et vecteur rouges).

Ainsi, la connaissance de l’ensemble des vecteurs vitesse du plan de phase, qu’on appelle le champ de vecteurs, permet d’obtenir l’ensemble des trajectoires solutions du système (Fig. 7).

Fig. 7 : Champ de vecteurs du système défini par \(\frac{dN}{dt} = f(N,S)\) et \(\frac{dS}{dt} = g(N,S)\).

Parmi tous les vecteurs vitesse du champ de vecteurs, certains sont remarquables, comme par exemple les vecteurs vitesse horizontaux et les vecteurs vitesse verticaux.

On définit :

• Les isoclines nulles horizontales, comme les courbes du plan de phase le long desquelles les vecteurs vitesse sont horizontaux. Un vecteur vitesse horizontal a pour coordonnées : \[\textrm{V} = \left( {\begin{array}{*{20}{c}} {\frac{dN}{dt}}\\ {0} \end{array}} \right)\] Par conséquent, les isoclines nulles horizontales sont les courbes solution de l’équation : \[\frac{dS}{dt} = 0 \Leftrightarrow g(N,S)=0\]

• Les isoclines nulles verticales, comme les courbes du plan de phase le long desquelles les vecteurs vitesse sont verticaux. Un vecteur vitesse vertical a pour coordonnées : \[\textrm{V} = \left( {\begin{array}{*{20}{c}} {0}\\ {\frac{dS}{dt}} \end{array}} \right)\] Par conséquent, les isoclines nulles horizontales sont les courbes solution de l’équation : \[\frac{dN}{dt} = 0 \Leftrightarrow f(N,S)=0\]

Ainsi, on peut compléter le champ de vecteurs avec les isoclines nulles dans le plan de phase (Fig. 8).

Fig. 8 : Représentation des isoclines nulles horizontales (en bleu) et verticales (en vert) du système défini par \(\frac{dN}{dt} = f(N,S)\) et \(\frac{dS}{dt} = g(N,S)\).

Puisque les points d’équilibre vérifient simultanément \(\frac{dN}{dt} = 0\) et \(\frac{dS}{dt} = 0\), on les trouve donc à l’intersection des isoclines horizontales et verticales (Fig. 8, points noirs).

3.6.2 Les trajectoires

Il ne reste plus maintenant qu’à dessiner quelques trajectoires pour différentes conditions initiales (Fig. 9).

library(phaseR)
m <- function(t,y,parameters)
 {   
  dy <- numeric(2)
  dy[1] <- parameters[1]*(y[2]/(1+y[2]))*y[1] - y[1]
  dy[2] <- -(y[2]/(1+y[2]))*y[1] - y[2] + parameters[2]
  list(dy)
}
# Points d'équilibre (expression analytique)
N2 <- a1*(a2-1/(a1-1))
S2 <- 1/(a1-1)
# Champ de vecteurs
par(mar=c(4, 4, 0.2, 0.2))
VectField <- flowField(m, xlim=c(0, 2.5), ylim=c(0, 2.5),
                       parameters=c(a1, a2), points=10, add=FALSE,
                       xlab="N(t)", ylab="S(t)", las=1, col="red")
grid()
# Isoclines
isocline <- nullclines(m, xlim=c(-0.1, 2.5), ylim=c(-0.1, 2.5),
                       parameters=c(a1, a2), points=150,
                       col=c("darkgreen","blue"), add=TRUE, lwd=2,
                       add.legend = FALSE)
# Points d'équilibre (sur le graphe)
points(c(0,N2), c(a2,S2), pch=19)
# Choix de conditions initiales
init1 <- c(0.5, 2.5)
trajectoire1 <- trajectory(m, y0=init1, tlim=c(0,100),
                           parameters=c(a1, a2), add=TRUE)
init2 <- c(1.5, 2.5)
trajectoire2 <- trajectory(m, y0=init2, tlim=c(0,100),
                           parameters=c(a1, a2), add=TRUE)
init3 <- c(0.5, 0)
trajectoire3 <- trajectory(m, y0=init3, tlim=c(0,100),
                           parameters=c(a1, a2), add=TRUE)
init4 <- c(1.5, 0)
trajectoire4 <- trajectory(m, y0=init4, tlim=c(0,100),
                           parameters=c(a1, a2), add=TRUE)
init5 <- c(2.5, 0)
trajectoire5 <- trajectory(m, y0=init5, tlim=c(0,100),
                           parameters=c(a1, a2), add=TRUE)

Fig. 9 : Trajectoires du système du système défini par \(\frac{dN}{dt} = f(N,S)\) et \(\frac{dS}{dt} = g(N,S)\).

3.6.3 Les chroniques

Comme évoqué précédemment, on peut aussi représenter les solutions d’un système dynamique dans les plans \((t,N(t))\) et \((t,S(t))\). On appelle ces représentations des chroniques, en référence au fait que ce sont des représentations au cours du temps (du grec,chronos qui signifie le temps) (Fig. 10).

par(mfrow=c(1,1), mar=c(4,4,0,0.5))
init1 <- c(0.5, 2.5)
par(mar=c(4, 4, 0.2, 0.2))
chronique <- numericalSolution(m, y0=init1, tlim=c(0,30),
                               parameters=c(a1, a2), type="one",
                               col=c("orange","gray"), las=1, lwd=2,
                               ylim=c(0,2.5), add.legend=FALSE,
                               xlab=expression(t),ylab="")
legend("topright", legend=c(expression(N(t)), expression(S(t))),
       col=c("orange","gray"), lwd=2, bty="n")

Fig. 10 : Chroniques associées au système défini par par \(\frac{dN}{dt} = f(N,S)\) et \(\frac{dS}{dt} = g(N,S)\).

3.6.4 Stabilité

Le logiciel R, et en particulier le package phaseR, permet de vérifier la stabilité des points d’équilibr.

• A l’aide du code R ci-dessous, testez différentes valeurs de paramètres. Justifiez d’un choix que vous conserverez pour la suite de cette partie.

stabilite <- stability(m, ystar = c(N2,S2), parameters= c(a1,a2),
                       system="two.dim", summary=FALSE)
stabilite$tr # critère stabilité 1

[1] -1.5

stabilite$Delta # critère stabilité 2

[1] 0.5

3.6.5 Points d’équilibre

Une fois le portrait de phase tracé, si on ne connaît pas l’expression analytique des points d’équilibre, il est possible de demander à R de les trouver par approximation numérique grâce à la fonction findEquilibrium() ; cela nécessite par contre d’avoir une idée approximative de l’endroit où ils se trouvent sur le portrait de phase, ce qui est possible si on a dessiné les isoclines.

eq1 <- findEquilibrium(m, y0=c(0.05,0.05), parameters=c(a1,a2),
                       summary=FALSE)
t(eq1$ystar) # Coordonnées du point d'équilibre

              [,1] [,2]
[1,] -1.440918e-09    2

# Stabilité du point d'équilibre ?
eq1$Delta # Critère 1

[1] -0.3333333

eq1$tr # Critère 2

[1] -0.6666667

• Faire de même pour l’autre point d’équilibre.

3.6.6 Analyse complète

• Explorer la fonction phasePlaneAnalysis().

3.7 Retour au modèle du chémostat

3.7.1 Écriture à 2 paramètres

Revenons maintenant au modèle du chémostat pour construire le portrait de phase et dessiner solutions et chroniques. On rappelle que le modèle s’écrit : \[\left\{ {\begin{array}{*{20}{l}} {\frac{{d{N}}}{{d{t}}} = {\alpha _1}\left( {\frac{{{S}}}{{1 + {S}}}} \right){N} - {N}}\\ {\frac{{d{S}}}{{d{t}}} = - \left( {\frac{{{S}}}{{1 + {S}}}} \right){N} - {S} + {\alpha _2}} \end{array}} \right.\]

Il faut en premier lieu choisir des valeurs de paramètres pour \(\alpha_1\) et \(\alpha_2\), sans oublier de respecter les deux conditions : \(\alpha_1 > 1\) et \(\alpha_2(\alpha_1 - 1) > 1\).

• Dessinez sous R le portrait de phase de ce système en ayant au préalable choisi un jeu de paramètres \((\alpha_1, \alpha_2)\) cohérents avec les figures présentées plus haut.

3.7.2 Retour au système d’équations initial

Revenons maintenant au modèle initial du chémostat avant la reparamétrisation (§ 3.4) : \[\left\{ \begin{array}{l} \frac{{dN\left( t \right)}}{{dt}} = \left( {\frac{{{r_{\max }}S\left( t \right)}}{{{K_s} + S\left( t \right)}}} \right)N\left( t \right) - \frac{F}{V}N\left( t \right)\\ \frac{{dS\left( t \right)}}{{dt}} = - \alpha \left( {\frac{{{r_{\max }}S\left( t \right)}}{{{K_s} + S\left( t \right)}}} \right)N\left( t \right) - \frac{F}{V}S\left( t \right) + \frac{F}{V}{S_0} \end{array} \right.\]

• Dessinez sous R le portrait de phase de ce système en ayant au préalable choisi un jeu de 6 paramètres cohérents expérimentalement.

Dans cette optique, vous utiliserez comme vecteur de paramètres pour encoder le modèle sous R, le vecteur suivant :

# parameters[1] = rmax
# parameters[2] = Ks
# parameters[3] = F
# parameters[4] = V
# parameters[5] = alpha
# parameters[6] = S0

• Retrouvez les coordonnées du point d’équilibre à partir des équations ci-dessus.

• Discutez de la signification des paramètres et de leur influence sur le comportement du système. Dans cet objectif, faites varier indépendamment chaque paramètre et regardez son influence sur le portrait de phase.

Eléments de correction pour le modèle du chémostat.