[INSA P2i-1] Bioréacteur & Modélisation

S. Charles, C. Lopes, J.R. Lobry & H. Charles

Contact : [Link]@[Link]

20 mai 2021
2
Contents

Accueil 5

1 Introduction 7
1.1 Généralités . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Quelques mots sur la modélisation . . . . . . . . . . . . . . . . . 9
1.3 Ce que nous allons faire . . . . . . . . . . . . . . . . . . . . . . . 10

2 Modélisation de la croissance des micro-organismes 13

2.1 Commençons en images . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Le modèle de croissance exponentielle (Malthus, 1798) . . . . . . 16
2.3 Le modèle de Verhulst (1845) . . . . . . . . . . . . . . . . . . . . 20
2.4 Introduction générale à R via RStudio . . . . . . . . . . . . . . . 25
2.5 Retour au modèle de Verhulst . . . . . . . . . . . . . . . . . . . . 26
2.6 Vers d’autres modèles . . . . . . . . . . . . . . . . . . . . . . . . 29
2.7 Comparaison des modèles de croissance bactérienne . . . . . . . . 32
2.8 Calcul du temps de doublement . . . . . . . . . . . . . . . . . . . 34

3 Modélisation du chémostat 35
3.1 Variables et paramètres . . . . . . . . . . . . . . . . . . . . . . . 36
3.2 Hypothèses et premières équations . . . . . . . . . . . . . . . . . 36
3.3 La notion de substrat limitant . . . . . . . . . . . . . . . . . . . . 37
3.4 Reparamétrisation du modèle . . . . . . . . . . . . . . . . . . . . 39
3.5 Analyse mathématique du modèle du chémostat . . . . . . . . . . 40
3.6 Représentations graphiques des solutions d’un système dynamique 43
3.7 Retour au modèle du chémostat . . . . . . . . . . . . . . . . . . . 50

4 Vitesse d’une réaction enzymatique 53

4.1 Les observations . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.2 La méthode de régression . . . . . . . . . . . . . . . . . . . . . . 54
4.3 Transformations linéaires . . . . . . . . . . . . . . . . . . . . . . 57
4.4 Retour à la régression gaussienne . . . . . . . . . . . . . . . . . . 60
4.5 A vous de jouer ! . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

3
4 CONTENTS
Accueil

Les auteurs de ce site :

• Sandrine CHARLES Page personnelle
• Christelle LOPES Page personnelle
• Jean LOBRY Page personnelle
• Hubert CHARLES Page personnelle

5
6 CONTENTS
Chapter 1

Introduction

1.1 Généralités

Un bioréacteur (Fig. 1) est une cuve dans laquelle se multiplient des micro-
organismes (levures, bactéries, champignons, algues…) qui consomment des sub-
strats ou se nourrissent d’autres organismes pour se développer, et qui utilisent
des précurseurs∗ et des activateurs∗ pour produire de la biomasse, synthétiser
des métabolites ou encore bioconvertir des molécules d’intérêt (e.g., dépollu-
tion). Grâce au bioréacteur, il est possible de contrôler les conditions de culture
(température, pH, aération…) et, de ce fait, de récolter des données expérimen-
tales relativement fiables pour le suivi de la croissance bactérienne et/ou de la
réaction chimique d’intérêt.

∗
Précurseur : composé chimique qui est consommé par une réaction (ou une
voie métabolique) et qui est transformé en un ou plusieurs autres composés.

∗
Activateur : substance qui, mélangée en proportion très faible, accélère la réac-
tion chimique (e.g., enzymes elles-mêmes ou co-facteurs enzymatiques comme
des ions ou des vitamines).

7
8 CHAPTER 1. INTRODUCTION

Fig. 1 : Schéma de principe d’un bioréacteur

On distingue trois modes d’alimentation classiques des bioréacteurs (Fig. 2) :
• Le mode “Batch” : à l’instant initial, la cuve est remplie par le milieu
de culture stérilisé (contenant substrats, précurseurs et activateurs), et les
espèces sont introduites. La croissance des micro-organismes se déroule
ensuite sans addition supplémentaire de milieu. Le volume reste constant.
On parle alors de fermenteur. En général, la productivité est faible. En
fin de processus, le fermenteur est vidé et son contenu est remplacé.
• Le mode “semi-continu” (ou “Fedbatch”) : à l’instant initial, la cuve est
remplie d’un volume initial 𝑉0 puis alimentée en continu par le milieu
de culture selon un débit réglé de façon à ce que la concentration en
substrat soit constante et ce jusqu’à son volume final. L’alimentation est
alors coupée. Le fedbatch permet en pratique un meilleur contrôle des
conditions de croissance, un gain de temps et la possibilité de modifier
le milieu en cours de culture. Il y a par contre un risque important de
contamination externe.
• Le mode “continu” (ou chémostat) : c’est le mode le plus largement
employé dans le domaine de la dépollution biologique de l’eau (Bernard,
2004). Caractérisé par un volume constant, un chémostat est soumis à un
1.2. QUELQUES MOTS SUR LA MODÉLISATION 9

soutirage du milieu de culture égal au taux d’alimentation de la cuve. Son

fonctionnement a été décrit dès 1950 par Novick et Szilard.

Fig. 2 : Les différents modes de fonctionnement d’un bioréeacteur construit par

l’homme ou d’un bioréacteur naturel (lac). D’après Bernard (2004).
Dans ce qui suit, nous considérerons des bioréacteurs infiniment mélangés (au
contenu supposé homogène) dont on peut modéliser le fonctionnement au moyen
d’équations différentielles ordinaires, c’est-à-dire en considérant le temps comme
une variable continue. Les variables modélisées sont généralement la biomasse
bactérienne (désignée par 𝑁 (𝑡) ci-après) ou la vitesse de la réaction enzymatique
(désignée par 𝑣 ci-après).

1.2 Quelques mots sur la modélisation

La modélisation des systèmes biologiques est une tâche extrêmement déli-
cate car, contrairement à la physique, il n’existe pas de lois admises et reconnues
par tous caractérisant l’évolution des phénomènes auxquels on s’intéresse.
La modélisation est la démarche scientifique qui permet l’élaboration d’un mod-
èle (Pavé, 2012). Elle est le plus souvent fondée sur les mathématiques, et même
la partie des mathématiques qui traite des variables et paramètres à valeurs
dans ℝ. Pourtant, d’autres formalismes peuvent parfois s’avérer utiles comme
les représentations schématiques dans les modèles de décision.
Modéliser, ce n’est pas théoriser. On peut modéliser sans théoriser et inverse-
ment. Par contre, un modèle est très souvent un outil précieux dans une dé-
marche théorique. De ce fait, la modélisation intervient dans les grandes fonc-
10 CHAPTER 1. INTRODUCTION

tions de la recherche scientifique :

• Détecter et énoncer des questions ;
• Transformer en problème et acquérir des données et des connaissances ;
• Définir des actions et étudier leurs conséquences.
Le modélisateur est celui qui modélise ; il est spécialiste d’une stratégie de
construction et d’utilisation de modèles : il maîtrise une grande variété de tech-
niques et de méthodes; il s’inspire du problème biologique pour proposer un
modèle (et non l’inverse) ; il s’implique dans la connaissance des aspects bi-
ologiques de ce qu’il modélise.
Le modélisateur modélise en proposant des modèles. Un modèle est une
représentation symbolique de certains aspects d’un phénomène du monde réel.
Ce n’est pas une fin en soi mais un outil parmi ceux de la boîte à outil du
modélisateur. Il est fortement couplé à l’expérience et/ou à l’observation. En
aucun cas il ne doit être le prétexte de décisions prises a priori ; un modèle ne
peut être qu’un instrument d’aide à la décision (technique ou politique). Par
définition, il sera donc toujours faux.
Pourtant, pour être efficace un modèle doit avant tout être opératoire, c’est-à-
dire permettre de répondre aux objectifs initiaux, être interprétable en termes
biologiques et être traduisible en termes simples et accessibles à tous.
L’élaboration d’un modèle doit donc prendre en compte :
• Le phénomène biologique à étudier ;
• Le formalisme choisi ;
• Les objectifs (que veut-on faire du modèle ?) ;
• Les données et les connaissances a priori, disponibles ou accessibles par
l’expérience ou l’observation.
Une fois le modèle écrit, le travail du modélisateur consiste à :
• Manipuler le modèle pour étudier ses propriétés ;
• Établir des relations avec d’autres représentations (graphes, simulations)
;
• Interpréter et confronter les résultats avec la réalité biologique, le plus
souvent vue au travers des données expérimentales.

1.3 Ce que nous allons faire

L’objectif de ce cours est donc de vous initier à la modélisation au travers de la
description des processus qui caractérisent le fonctionnement d’un bioréacteur.
Nous nous intéresserons d’abord à la croissance des micro-organismes présents
1.3. CE QUE NOUS ALLONS FAIRE 11

à l’intérieur d’un bioréacteur ; puis nous examinerons la dynamique couplée de

ces micro-organismes à celle de la quantité de substrat disponible ; enfin, nous
nous intéresserons à la vitesse d’une réaction enzymatique à l’image de celles
que l’on trouve au sein d’un bioréacteur.
Ces mises en situation biologiques nous permettrons d’aborder des aspects plus
méthodologiques : la résolution exacte d’équations différentielles ordinaires,
l’étude qualitative de systèmes dynamiques et la régression non linéaire, c’est-
à-dire l’ajustement d’un modèle à des données expérimentales pour en estimer
les paramètres. Ces points seront traités conjointement par l’utilisation d’outils
mathématiques et de simulations numériques qui seront réalisées à l’aide du
logiciel R et de son interface utilisateur RStudio.
Pour éviter au maximum de passer du temps sur les aspects techniques, nous
vous proposons ce cours sous la forme d’un tutoriel à suivre pas à pas ; il intègre
directement les outils R dont vous aurez besoin sous la forme d’une interface
Shiny.
12 CHAPTER 1. INTRODUCTION
Chapter 2

Modélisation de la
croissance des
micro-organismes

La notion de croissance bactérienne recouvre deux aspects : (1) la croissance

de la cellule bactérienne elle-même (taille, poids, volume) ; (2) le phénomène de
division cellulaire (croissance de la population bactérienne) auquel nous nous
intéresserons exclusivement ici.

2.1 Commençons en images

2.1.1 La croissance bactérienne
• La fission binaire

Les bactéries sont des êtres vivants unicellulaires qui se reproduisent par fission
binaire (ou scissiparité) : une cellule mère donne deux cellules filles. Ces deux
cellules filles donnent ensuite quatre cellules petite-filles…

Pour en savoir plus…(consulté le 2021-05-01)

• La croissance bactérienne

13
14CHAPTER 2. MODÉLISATION DE LA CROISSANCE DES MICRO-ORGANISMES

La croissance de la population bactérienne correspond à la répétition du proces-

sus de reproduction bactérienne décrit ci-dessus.

• Mesure de la croissance bactérienne

La croissance bactérienne peut se mesurer par deux méthodes :

• soit par une mesure de Densité Optique (DO) ; au laboratoire, cette mesure
de turbidimétrie est la plus pratique et rapide. On choisit en général une
longueur d’onde dans le visible, soit 600 nm. La densité optique de la
solution croît linéairement avec le nombre de bactéries.

• soit par un dénombrement cellulaire sur boîte de Pétri (voir figure ci-
dessous).

D’après ce document (consulté le 2021-05-01).

Lors d’un dénombrement sur boîte de Petri, on prend en considération les boites
de 30 à 300 colonies. Le calcul se fait comme suit : si pour la dilution 1/10000,
on obtient 125 colonies, alors on dira que l’on a 125 x 10000= 1,25 million de
bactéries par ml. On suppose dans ce cas que chaque colonie est issue d’une
bactérie, mais en réalité c’est faux. On parle plutôt d’Unités Formant Colonies
(UFC), chaque colonie provenant de plusieurs bactéries en chainettes ou en amas.
2.1. COMMENÇONS EN IMAGES 15

2.1.2 Les différentes phases de la croissance bactérienne

C’est Buchanan, en 1918, qui identifie le premier que la croissance bactérienne
se décompose en plusieurs phases. Monod, en 1949, caractérisera les six phases
en fonction du taux de croissance, comme le montre la figure ci-dessous (Fig.
2bis) :
1. La phase de latence, pendant laquelle le taux de croissance est nul ;
2. La phase d’accélération, pendant laquelle le taux de croissance augmente
;
3. La phase exponentielle, pendant laquelle le taux de croissance est constant
et maximum ;
4. La phase de ralentissement, pendant laquelle le taux de croissance diminue
;
5. La phase stationnaire, pendant laquelle le taux de croissance est redevenu
nul ;
6. La phase de déclin, pendant laquelle le taux de croissance est négatif.

Fig. 2bis : les différentes phase de la croissance bactérienne (Monod, 1949)

Dans ce cours nous nous focaliserons uniquement sur les phase 2 à 5.
16CHAPTER 2. MODÉLISATION DE LA CROISSANCE DES MICRO-ORGANISMES

2.1.3 La dynamique des populations

Images mathématiques de la Lutte pour la Vie Images mathématiques de la
Lutte pour la Vie est un film de Jean Painlevé (1937) qui évoque les lois de la
croissance et de l’hérédité des populations (visionnage optionnel, durée 13 min
30).

2.2 Le modèle de croissance exponentielle

(Malthus, 1798)

2.2.1 Modélisation de la dynamique d’une population

L’accroissement d’une population isolée est en première approximation fonction
des naissances, des morts et des processus de migration des individus de la
population vers ou depuis un autre environnement. Bien que suggéré très tôt
par Euler (≃ 1760), on attribue cependant à Sir Thomas Robert Malthus (1798)
le modèle le plus simple proposé pour décrire l’évolution dans le temps d’une
population isolée.
Les hypothèses sous-jacentes au modèle de Malthus sont les suivantes :
• Le processus de migration est négligé (système fermé) ;
• L’accroissement absolu de la population en termes d’effectif (resp. den-
sité ou biomasse) est supposé proportionnel à l’effectif (resp. densité ou
biomasse), et à la longueur de l’intervalle de temps, selon une échelle con-
tinue, pendant lequel on mesure cet accroissement ;
• Les individus de la population sont supposés isolés ou équivalents, i.e.,
qu’on ne prend pas en compte ni d’interaction entre individus, ni de struc-
ture d’âge, ni d’auto-régulation de la croissance ;
2.2. LE MODÈLE DE CROISSANCE EXPONENTIELLE (MALTHUS, 1798)17

• La taille de la population (en termes d’effectif, de densité ou de biomasse)

est correctement représentée par sa moyenne.

La traduction directe de ces hypothèses conduit à écrire :

Δ𝑁 (𝑡) = 𝑟𝑁 (𝑡)Δ𝑡

où 𝑁 (𝑡) représente l’effectif (resp. la densité ou la biomasse) de la population au

temps 𝑡, Δ𝑁 (𝑡) l’accroissement absolu de la population, Δ𝑡 l’intervalle de temps
pendent lequel on mesure l’accroissement, et 𝑟 le coefficient de proportionnalité.

On considérera que 𝑁 (𝑡) > 0.

En supposant que le raisonnement reste valable pour de petites variations de 𝑡

(Δ𝑡 → 0), il vient :

𝑑𝑁 (𝑡)
= 𝑟𝑁 (𝑡)
𝑑𝑡

1 𝑑𝑁(𝑡)
Ainsi, le coefficient 𝑟 = 𝑁(𝑡) 𝑑𝑡 représente le taux de croissance relatif (ou
intrinsèque) de la population, que l’on appelle encore taux de croissance
malthusien. Ce paramètre 𝑟 intègre donc à la fois les naissances et les morts
qui influencent la dynamique de la population :

𝑟 =𝑏−𝑑

avec 𝑏 > 0 le taux de natalité naturelle de la population et 𝑑 > 0 le taux de

mortalité naturelle.

En première intention le modèle de Malthus fût proposé pour décrire l’évolution

de la population mondiale (Fig. 3) :
18CHAPTER 2. MODÉLISATION DE LA CROISSANCE DES MICRO-ORGANISMES

5
Population mondiale

1750 1800 1850 1900 1950 2000

Année

Fig. 3 : Evolution de la population mondiale (exprimée en milliards d’habitants)

de 1750 à nos jours (données manquantes entre 1750 et 1900).
2.2. LE MODÈLE DE CROISSANCE EXPONENTIELLE (MALTHUS, 1798)19

2.2.2 Retour à la croissance bactérienne

Laboratoire virtuel

2.2.3 Analyse du modèle de Malthus

Le modèle de Malthus est une équation différentielle ordinaire, d’ordre 1 et à
variables séparables (ici 𝑡 et 𝑁 (𝑡)), qui peut donc s’intégrer facilement :

𝑑𝑁(𝑡)
𝑑𝑡 = 𝑟𝑁 (𝑡)
⇔ 𝑑𝑁(𝑡)
𝑁(𝑡) = 𝑟𝑑𝑡
⇔ ∫ 𝑑𝑁(𝑡)
𝑁(𝑡) = ∫ 𝑟𝑑𝑡
⇔ ln 𝑁 (𝑡) = 𝑟𝑡 + 𝑐𝑠𝑡𝑒
⇔ 𝑁 (𝑡) = 𝐶𝑒𝑟𝑡

Finalement 𝑁 (𝑡) = 𝑁0 𝑒𝑟𝑡 avec 𝑁0 l’inocculum bactérien initial. La solution

du modèle de Malthus est une fonction exponentielle, d’où son nom de modèle
exponentiel. De plus :

• Si 𝑟 > 0 (i.e., 𝑏 > 𝑑), alors lim 𝑁 (𝑡) = +∞ ; le taux de natalité est plus
𝑡→+∞
élevé que le taux de mortalité, la population croît à l’infini.

• Si 𝑟 = 0 (i.e., 𝑏 = 𝑑), alors ∀𝑡 𝑁 (𝑡) = 𝑁0 ; le taux de natalité est égal au

taux de mortalité, la population reste stable.

• Si 𝑟 < 0 (i.e., 𝑏 < 𝑑), alors lim 𝑁 (𝑡) = 0 ; le taux de natalité est moins
𝑡→+∞
élevé que le taux de mortalité, la population s’éteint.

• D’après vous, quels sont les limitations de ce modèle ? Revenir sur les
hypothèses de départ.
20CHAPTER 2. MODÉLISATION DE LA CROISSANCE DES MICRO-ORGANISMES

2.3 Le modèle de Verhulst (1845)

2.3.1 L’équation logistique

Même si le modèle de Malthus est encore utilisé pour décrire certains
phénomènes biologiques, il n’en reste pas moins biologiquement irréaliste par le
fait qu’il repose sur l’hypothèse d’une croissance sans aucune limitation ni par
l’espace, ni par la ressource, ni par la densité d’individus dans l’environnement
considéré. Des modèles alternatifs ont donc été proposés, introduisant en
particulier la notion d’autorégulation de la population, en tenant compte
des phénomènes de densité-dépendance ; on parle aussi de compétition
intra-spécifique.

C’est ainsi que Pierre-François Verhulst, mathématicien belge, proposa en 1845

son fameux modèle logistique. Les hypothèses sous-jacentes à ce modèle sont
les suivantes :

• Le taux de natalité dépend de la taille de la population selon une relation

linéaire décroissante : 𝑏(𝑁 ) = 𝑏0 − 𝛽𝑁 (𝑡), avec 𝑏0 le taux de natalité
naturelle et 𝛽 un coefficient de densité-dépendance ; 𝑏0 , 𝛽 > 0.

• Le taux de mortalité dépend de la taille de la population selon une relation

linéaire croissante : 𝑑(𝑁 ) = 𝑑0 + 𝛿𝑁 (𝑡), avec 𝑑0 le taux de mortalité
naturelle et 𝛿 un coefficient de densité-dépendance ; 𝑑0 , 𝛿 > 0.

Le modèle de Verhulst suppose de plus que 𝑏0 > 𝑑0 , i.e. lorsque 𝑁 est petit, il
y a effectivement croissance de la population.
2.3. LE MODÈLE DE VERHULST (1845) 21

Si on repart du modèle de Malthus, il vient :

𝑑𝑁(𝑡)
𝑑𝑡 = (𝑏 (𝑁 (𝑡)) − 𝑑 (𝑁 (𝑡))) 𝑁 (𝑡)
⇔ 𝑑𝑁(𝑡)
𝑑𝑡 = (𝑏0 − 𝛽𝑁 (𝑡) − (𝑑0 + 𝛿𝑁 (𝑡))) 𝑁 (𝑡)
𝑑𝑁(𝑡)
⇔ 𝑑𝑡 = (𝑏0 − 𝑑0 ) 𝑁 (𝑡) − (𝛽 + 𝛿) 𝑁 2 (𝑡)
⇔ 𝑑𝑁(𝑡) 𝛽+𝛿
𝑑𝑡 = (𝑏0 − 𝑑0 ) 𝑁 (𝑡) (1 − 𝑏0 −𝑑0 𝑁 (𝑡))
⇔ 𝑑𝑁(𝑡) 𝑁(𝑡)
𝑑𝑡 = 𝑟0 𝑁 (𝑡) (1 − 𝐾 )

où 𝑟0 = 𝑏0 − 𝑑0 est le taux de croissance malthusien de la population (𝑟0 > 0)

et 𝐾 = 𝑏0𝛽+𝛿
−𝑑0
désigne la capacité limite de la population (𝐾 > 0). Nous
reviendrons plus tard sur cette notion.
Cette équation est le fondement du modèle évolutif 𝑟/𝐾, c’est-à-dire la clas-
sification des espèces en fonction de leur stratégie de reproduction : stratégie
𝑟 pour les espèces qui produisent beaucoup de jeunes avec une forte mortalité
juvénile (𝑟 grand, 𝐾 petit); stratégie 𝐾 pour les espèces qui produisent peu de
jeunes avec une faible mortalité juvénile (𝑟 faible, 𝐾 grand). Cette équation
sera étendue au cas de deux populations en compétition un siècle plus tard par
le mathématicien italien Vito Volterra.

2.3.2 Retour à la croissance bactérienne

Reprenons les expériences virtuelles de croissance bactérienne mais sur une péri-
ode de temps cette fois beaucoup plus longue.
Laboratoire virtuel

2.3.3 Analyse du modèle de Verhulst

[Link] Analyse qualitative
Pour analyser le modèle de Verhulst, il est commode de décomposer la variation
dans le temps de la manière suivante :
𝑑𝑁 (𝑡) 𝑁 2 (𝑡)
= 𝑟0 𝑁 (𝑡) − 𝑟0
𝑑𝑡 𝐾
On comprend ainsi que la croissance naturelle de la population, représentée par
le terme 𝑟0 𝑁 (𝑡), est régulée négativement par le terme de compétition intra-
2
spécifique −𝑟0 𝑁𝐾(𝑡) (cf. film de J. Painlevé), qui prévaut d’autant plus que 𝑁 (𝑡)
est grand.
Le terme de compétition intra-spécifique (ou compétition de Verhulst) influ-
ence d’autant plus négativement la croissance naturelle de la population que le
paramètre 𝐾 est petit :
• Si 𝐾 → +∞, alors le phénomène de densité-dépendance est négligeable ;
• Si 𝐾 → 0, alors le phénomène de densité-dépendance est prépondérant et
contraint fortement la croissance de la population.
22CHAPTER 2. MODÉLISATION DE LA CROISSANCE DES MICRO-ORGANISMES

C’est pourquoi on appelle le paramètre 𝐾 la capacité limite de la population,

qui correspond en fait à une capacité limité d’accueil du milieu dans lequel
évolue la population.

Du point de vue dynamique, on imagine alors ce qu’il peut se passer :

• Si 𝑁 (𝑡) est petit devant 𝐾, alors 𝑑𝑁(𝑡)

𝑑𝑡 ≃ 𝑟0 𝑁 (𝑡). Comme 𝑟0 > 0, la
population va croître de manière exponentielle ;

• Lorsque 𝑁 (𝑡) ≃ 𝐾 (c’est-à-dire 𝑁 (𝑡) proche de 𝐾), il peut se passer deux

choses :

1. Si 𝑁 (𝑡) < 𝐾, alors on peut écrire :

𝑑(𝐾−𝑁(𝑡))
𝑑𝑡 = − 𝑑𝑁(𝑡)
𝑑𝑡
⇔ 𝑑(𝐾−𝑁(𝑡))
𝑑𝑡 = −𝑟 𝑁(𝑡)
0 𝑁 (𝑡) (1 − 𝐾 )
𝑑(𝐾−𝑁(𝑡)) 𝑁(𝑡)
⇔ 𝑑𝑡 = −𝑟0 𝐾 (𝐾 − 𝑁 (𝑡))

𝑁(𝑡)
avec 𝐾 ≃ 1.

On retrouve ainsi un modèle de type exponentiel pour la variable 𝐾 −𝑁 (𝑡) avec

un taux ≃ −𝑟0 puisque 𝑁 (𝑡) ≃ 𝐾 (i.e., 𝑁(𝑡)
𝐾 ≃ 1). Ce taux est négatif ce qui
veut dire que la limite de la variable 𝐾 − 𝑁 (𝑡) lorsque 𝑡 tend vers +∞ est 0,
c’est-à-dire 𝑁 (𝑡) → 𝐾.

2. Si 𝑁 (𝑡) > 𝐾, alors on peut écrire :

𝑑(𝑁(𝑡)−𝐾)
𝑑𝑡 = 𝑑𝑁(𝑡)
𝑑𝑡
⇔ 𝑑(𝑁(𝑡)−𝐾)
𝑑𝑡 = 𝑟0 𝑁 (𝑡) (1 − 𝑁(𝑡)
𝐾 )
𝑑(𝑁(𝑡)−𝐾) 𝑁(𝑡)
⇔ 𝑑𝑡 = −𝑟0 𝐾 (𝑁 (𝑡) − 𝐾)

𝑁(𝑡)
avec 𝐾 ≃ 1.

On retrouve comme précédemment un modèle de type exponentiel, mais cette

fois-ci pour la variable 𝑁 (𝑡) − 𝐾 dont le taux ≃ −𝑟0 est < 0, ce qui conduit de
nouveau à 𝑁 (𝑡) → 𝐾.

En conséquence, quelque soit la condition initiale 𝑁 (0) > 0, on aura toujours

lim 𝑁 (𝑡) = 𝐾. Il y a donc bien auto-régulation de la dynamique de la pop-
𝑡→+∞
ulation, dont l’effectif (respectivement la densité ou la biomasse) ne peut pas
dépasser la capacité limite 𝐾.
2.3. LE MODÈLE DE VERHULST (1845) 23

[Link] Résolution exacte

De nouveau, le modèle de Verhulst est une équation différentielle ordinaire
d’ordre 1 à variables séparables que l’on sait résoudre explicitement :
𝑑𝑁(𝑡) 𝑁(𝑡)
𝑑𝑡 = 𝑟0 𝑁 (𝑡) (1 − 𝐾 )
𝑟0
⇔ 𝑑𝑁(𝑡)
𝑑𝑡 = 𝐾 𝑁 (𝑡) (𝐾 − 𝑁 (𝑡))
⇔ 𝑁(𝑡)(𝐾−𝑁(𝑡)) = 𝑟𝐾0 𝑑𝑡
𝑑𝑁(𝑡)

𝑑𝑁(𝑡)
⇔ ∫ 𝑁(𝑡)(𝐾−𝑁(𝑡)) = ∫ 𝑟𝐾0 𝑑𝑡
𝑑𝑁(𝑡)
⇔ ∫ 𝑁(𝑡)(𝐾−𝑁(𝑡)) = 𝑟𝐾0 𝑡 + 𝑐𝑠𝑡𝑒
Pour continuer et intégrer le terme de gauche de cette équation, il faut d’abord
faire une décomposition en éléments simples :
1 𝐴 𝐵
𝑁(𝑡)(𝐾−𝑁(𝑡)) = 𝑁(𝑡) + 𝐾−𝑁(𝑡)
1
⇔ 𝑁(𝑡)(𝐾−𝑁(𝑡)) = 𝐴(𝐾−𝑁(𝑡))+𝐵𝑁(𝑡)
𝑁(𝑡)(𝐾−𝑁(𝑡))
1
⇔ 𝑁(𝑡)(𝐾−𝑁(𝑡)) = 𝐴𝐾+(𝐵−𝐴)𝑁(𝑡)
𝑁(𝑡)(𝐾−𝑁(𝑡))

ce qui implique :
𝐴𝐾 = 1
{
𝐵−𝐴=0
𝐴 = 1/𝐾
⇔{
𝐵=𝐴

et donc :
1 1 1 1
= ( + )
𝑁 (𝑡) (𝐾 − 𝑁 (𝑡)) 𝐾 𝑁 (𝑡) 𝐾 − 𝑁 (𝑡)

Ainsi :
𝑑𝑁(𝑡)
∫ 𝑁(𝑡)(𝐾−𝑁(𝑡)) =𝐾1
(∫ 𝑑𝑁(𝑡) 𝑑𝑁(𝑡)
𝑁(𝑡) + ∫ 𝐾−𝑁(𝑡) )
1
= 𝐾 (ln (𝑁 (𝑡)) − ln (𝐾 − 𝑁 (𝑡)))
1 𝑁(𝑡)
=𝐾 ln ( 𝐾−𝑁(𝑡) )

Finalement, on obtient :
1 𝑁(𝑡) 𝑟0
𝐾 ln ( 𝐾−𝑁(𝑡) )= 𝐾𝑡 + 𝑐𝑠𝑡𝑒
𝑁(𝑡)
⇔ ln ( 𝐾−𝑁(𝑡) ) = 𝑟0 𝑡 + 𝑐𝑠𝑡𝑒
𝑁(𝑡) 𝑟0 𝑡
⇔ 𝐾−𝑁(𝑡) = 𝐶𝑒
𝑟0 𝑡
⇔ 𝑁 (𝑡) = 𝐶𝑒 (𝐾 − 𝑁 (𝑡))
⇔ 𝑁 (𝑡) (1 + 𝐶𝑒𝑟0 𝑡 ) = 𝐶𝐾𝑒𝑟0 𝑡
𝐶𝐾𝑒𝑟0 𝑡
⇔ 𝑁 (𝑡) = 1+𝐶𝑒 𝑟0 𝑡

Il reste à déterminer la constante d’intégration 𝐶 à partir de la condition initiale

𝑁 (0) = 𝑁0 , ce qui conduit à :
𝐶𝐾 𝑁0
𝑁0 = ⇔ 𝐶𝐾 = 𝑁0 (1 + 𝐶) ⇔ 𝐶 =
1+𝐶 𝐾 − 𝑁0
24CHAPTER 2. MODÉLISATION DE LA CROISSANCE DES MICRO-ORGANISMES

La solution exacte du modèle de Verhulst s’écrit donc finalement :

𝑁0 )𝐾𝑒𝑟0 𝑡
( 𝐾−𝑁
𝑁 (𝑡) = 0
𝑁
0 )𝑒𝑟0 𝑡
1+( 𝐾−𝑁
0
𝑁0 𝐾𝑒𝑟0 𝑡
⇔𝑁 (𝑡) = (𝐾−𝑁 𝑟0 𝑡
0 )+𝑁0 𝑒
⇔𝑁 (𝑡) = (𝐾−𝑁 𝑁)𝑒0−𝑟𝐾
𝑡
0 +𝑁
0 0
𝐾
⇔𝑁 (𝑡) = 𝐾
1+( 𝑁 −1)𝑒−𝑟0 𝑡
0

On retrouve le résultat établi précédemment, à savoir que lim 𝑁 (𝑡) = 𝐾

𝑡→+∞
puisque lim 𝑒−𝑟0 𝑡 = 0.
𝑡→+∞

On peut par ailleurs démontrer que lorsque 𝑁 (𝑡) est petit et 𝑁0 << 𝐾, alors
𝑁 (𝑡) ≃ 𝑁0 𝑒𝑟0 𝑡 (*), c’est-à-dire que lorsque 𝑁 (𝑡) est petit, la croissance est
exponentielle avec un taux de croissance malthusien égale à 𝑟0 .

(*) Si 𝑁0 << 𝐾, alors 𝐾/𝑁0 >> 1 et (𝐾/𝑁0 − 1)𝑒−𝑟0 𝑡 ∼ (𝐾/𝑁0 )𝑒−𝑟0 𝑡

[Link] Résumer l’information : le portrait de phase

Nous avons démontré analytiquement à partir de la solution du modèle logistique
et vérifié par les simulations que lim 𝑁 (𝑡) = 𝐾.
𝑡→+∞

Revenons à l’équation différentielle ordinaire du modèle :

𝑑𝑁 (𝑡) 𝑁 (𝑡)
= 𝑟0 𝑁 (𝑡) (1 − )
𝑑𝑡 𝐾

Si l’on cherche les solutions constantes de cette équation, qu’on appelle des
points d’équilibre, c’est-à-dire les solutions 𝑁 (𝑡) telles que 𝑑𝑁(𝑡)
𝑑𝑡 = 0, il vient
𝑁1∗ = 0 et 𝑁2∗ = 𝐾.

En dehors de ces points d’équilibre, les solutions 𝑁 (𝑡) sont soient croissantes
( 𝑑𝑁(𝑡) 𝑑𝑁(𝑡)
𝑑𝑡 > 0) soit décroissantes ( 𝑑𝑡 < 0).

On peut donc résumer entièrement la dynamique du modèle logistique sur ce

que l’on appelle un portrait de phase :
2.4. INTRODUCTION GÉNÉRALE À R VIA RSTUDIO 25

2.4 Introduction générale à R via RStudio

• Ouvrir le logiciel RStudio.

RStudio (consulté le 2021-05-01) est un environnement de développement dédié

au logiciel libre de statistiques R (consulté le 2021-05-01). RStudio permet de
voir l’ensemble des objets créés dans son espace de travail. Il simplifie l’accès
à l’historique des commandes. Il facilite l’installation de nouveaux paquets
(consulté le 2021-05-01) et le chargement de ceux-ci dans l’espace de travail.
Il permet aussi de naviguer à travers tous les graphiques créés au cours de la
session de travail et de zoomer sur chacun de ces graphiques par un simple
clic. L’éditeur de programmes gère la coloration syntaxique, la complétion de
commandes et même la complétion de noms de fichiers (exactement comme
dans un terminal). Enfin, avec RStudio, chacun peut avoir accès à un tableur
permettant de parcourir ses données facilement (d’après framasoft (consulté le
2021-05-01)).

2.4.1 Quelques commandes de base avec R

getwd() # renvoie le chemin d'accès complet au dossier
# dans lequel vous vous trouvez

Vérifiez que c’est bien votre dossier de travail.

ls() # renvoie la liste des objets crés dans l'espace de travail

?nom-fonction # : renvoie la documentation concernant

# la fonction 'nom-fonction'
# Elle peut s'obtenir aussi avec la commande :
help("nom-fonction") # N'hésitez pas à en abuser

library() # renvoie la liste des librairies installées

# et donc disponibles.
library(nom-library) # charge la librairie spéficiée par 'nom-library'
library(help="nom-library")} # liste les fonctions disponibles
# au sein de la libraire 'nom-library'

'#' # permet d'insérer un commentaire

Voici quelques autres commandes utiles (pensez à visiter l’aide avec ?) :

mode, str, [Link], [Link] # types d'objets
length, dim, names # dimensions d'objets
plot, curve # fonctions graphiques de niveau 1
curve, lines, points # fonctions graphiques de niveau 2
26CHAPTER 2. MODÉLISATION DE LA CROISSANCE DES MICRO-ORGANISMES

2.4.2 Les objets R

== # comparaison entre les termes à droite et à gauche
# Par exemple
x <- 7 # affecte la valeur 7 à l'objet x
# De même : 7 -> x ou encore x = 7

Les valeurs manquantes sont représentées par ‘NA’. L’indétermination se

propage systématiquement dans les calculs.
Les vecteurs sont des objets au contenu homogène alors que les listes sont des
objets au contenu hétérogène. Les dataframe sont des listes dont tous les
éléments ont la même longueur (tableau).
x <- 1:5 # génère une série d'entiers de 1 à 5
y <- c("a","b","c") # génère une combinaison (i.e., un vecteur)
# des valeurs a, b et c
y[2] # renvoie le deuxième élément de y
s1 <- rep(1:5,2) # répète 2 fois la série d'entiers de 1 à 5
s2 <- seq(from=1, to=5.5, by=0.5) # génère une série de nombres
# allant de 1 à 5 avec un pas de 0.5
l <- list(s1,s2) # renvoie une liste à deux éléments s1 et s2
# On accède à un élément d'une liste par son nom à l'aide de
l[[1]]
d <- [Link](s1,s2) # renvoie un tableau à deux colonnes et 10 lignes
d
# On accède à un élément d'un tableau à l'aide des crochets :
d[1, 2]
d$s1 # ou d[,1] renvoie la première colonne du tableau d
d[1,] # renvoie la première ligne du tableau d

2.4.3 Pour aller plus loin

Prise en main de RStudio (consulté le 2021-05-01)
A brief introduction to R (consulté le 2021-05-01)

2.5 Retour au modèle de Verhulst

2.5.1 Simulation d’une chronique
Considérons le cas de la bactérie E. coli à 20°C.
• Copier-coller les lignes de code ci-dessous et observez ce qui se passe.
NB : Avec RStudio, créez un nouveau fichier de type ‘R script’ que vous enreg-
istrerez dans votre espace de travail (‘script_nom.R’). Incrémenterez ce fichier
au fur et à mesure de votre exploration des commandes R.
2.5. RETOUR AU MODÈLE DE VERHULST 27

N0 <- 100
K <- 9000
r0 <- 0.85
# courbe représentative d'une fonction
# x représente ici le temps
par(mar=c(4, 4, 0.2, 0.2))
curve(K/(1+(K/N0-1)*exp(-r0*x)), from=0, to=15, las=1,
main="Le modèle logistique", xlab="Temps (h)",
ylab="Taille de la population bactérienne (UFC)")

• A quoi correspondent les options from et to ?

• Que fait l’option las=1 ?

• Utilisez l’aide de la fonction curve pour explorer le comportement du mod-

èle logistique :

– En gardant 𝑁0 et 𝑟0 à leur valeur par défaut, faites varier les valeurs

de 𝐾 et superposer les différentes courbes en leur attribuant des
couleurs différentes.

Que remarquez-vous ?

– En gardant 𝑁0 et 𝐾 à leur valeur par défaut, faites varier les valeurs

de 𝑟0 et superposer les différentes courbes en leur attribuant des
couleurs différentes.

Que remarquez-vous ?

– En gardant 𝑟0 et 𝐾 à leur valeur par défaut, faites varier les valeurs de

𝑁0 (vous essaierez en particulier des valeurs 𝑁0 > 𝐾) et superposer
les différentes courbes en leur attribuant des couleurs différentes.

Que remarquez-vous ?

2.5.2 Comparaison avec les données observées

• En vous s’inspirant des instructions ci-dessous, représentez graphiquement
vos données sauvegardées

(Utilisez vos fichiers “[Link]” et/ou “[Link]”).

data <- [Link]("[Link]", header=TRUE, sep="\t")
par(mar=c(4, 4, 0.2, 0.2)) # réglage des marges
plot(data$Temps, data$Nbcol, xlim=c(0,5), ylim=c(0,max(data$N)),
pch=19, las=1)
28CHAPTER 2. MODÉLISATION DE LA CROISSANCE DES MICRO-ORGANISMES

6000

5000

4000
data$Nbcol

3000

2000

1000

0 1 2 3 4 5

data$Temps

• Peaufiner l’esthétique de votre graphique :

– Changez le type de point
– Labellisez les axes
– Mettez un titre au graphique
– Changez la couleur des points
– Ajoutez une légende
• Superposer un modèle à vos données expérimentales.
Vous choisirez le modèle exponentiel ou le modèle logistique selon les données
que vous avez représentées.
Vous choisirez un jeu de paramètres qui convient, en justifiant de votre choix.

2.5.3 Portrait de phase

Comme évoqué plus haut, on peut donc résumer entièrement la dynamique du
modèle logistique avec le portrait de phase, qui peut être construit avec R :
par(mar=c(0, 0, 0, 0))
plot(0,0, xlim=c(-1,1), ylim=c(-0.25,0.25), xaxt="n", yaxt="n",
xlab="", ylab="", bty="n")
text(1,0.1, labels=expression(N(t)))
arrows(-1,0,1,0)
points(c(-1,0), c(0,0), pch=19)
text(c(-1,0), c(-0.1,-0.1), labels=c("0","K"))
arrows(-0.6,0,-0.5,0, length=0.1, col="red", lwd=2)
2.6. VERS D’AUTRES MODÈLES 29

arrows(0.6,0,0.5,0, length=0.1, col="red", lwd=2)

par(mar=c(5,4,4,2)) # retour aux valeurs par défaut

Générez la figure correspondant au code R ci-dessus.

Revenons à la courbe représentative du modèle logistique. On remarque un
changement de courbure : pour des temps petits, la courbe est convexe, pour
les temps plus grands, la courbe est concave. Autrement dit, il existe un point
2
d’inflexion que l’on peut caractériser et qui vérifie 𝑑 𝑑𝑡
𝑁(𝑡)
2 = 0.
𝑑𝑁(𝑡) 𝑁(𝑡)
Posons 𝑑𝑡 = 𝑓 (𝑁 (𝑡)) avec 𝑓(𝑁 ) = 𝑟0 𝑁 (𝑡) (1 − 𝐾 ).

𝑑2 𝑁(𝑡) 𝑑 𝑑𝑓(𝑁(𝑡)) 𝑑𝑁(𝑡)

Ainsi 𝑑𝑡2 = 𝑑𝑡 (𝑓 (𝑁 (𝑡))) = 𝑑𝑁 × 𝑑𝑡 .

Par conséquent, le point d’inflexion n’étant pas un point d’équilibre (ce qui
impose 𝑑𝑁(𝑡)
𝑑𝑡 ≠ 0), il est solution de :

𝑑𝑓(𝑁(𝑡))
𝑑𝑁 =0
⇔ 𝑟0 (1 − 2𝑁(𝑡)
𝐾 )=0
⇔ 𝑁inf = 𝐾2

En reprenant votre code R précédent, symbolisez sur la courbe représentative

du modèle logistique, le point d’inflexion par une droite horizontale en pointillés.
Pour cela, vous utiliserez la fonction abline.
Eléments de correction pour le modèle de Verhulst.

2.6 Vers d’autres modèles

Il existe bien d’autres modèles concurrents du modèle logistique pour décrire
la croissance d’une population bactérienne (voir Pavé (2012), p126-127). Nous
vous proposons d’en explorer deux en particulier, un troisième est évoqué pour
information. Pour les étudier, vous suivrez, pour chacun des deux premiers, les
étapes suivantes :
• Résolution exacte pour la condition initiale 𝑁 (0) = 𝑁0 avec 0 < 𝑁0 < 𝐾
;
• Simulations pour différentes valeurs de 𝑟0 et 𝐾. On prendra 𝑁0 = 100 ;
• Recherche des points d’équilibre et de leur stabilité ;
• Détermination du point d’inflexion.

2.6.1 Le modèle de Gompertz

Dans l’équation proposée par Verhulst, le terme 1 − 𝑁 (𝑡)/𝐾 correspond au
freinage logistique évoqué dans le film de Jean Painlevé. C’est en effet ce
30CHAPTER 2. MODÉLISATION DE LA CROISSANCE DES MICRO-ORGANISMES

terme qui, multiplié au terme de croissance exponentiel, permet d’obtenir la

stabilisation de la croissance à la limite asymptotique 𝐾.
Benjamin Gompertz, en 1825, proposa un terme de freinage un peu différent
pour aboutir au modèle suivant :

𝑑𝑁 (𝑡) 𝐾
= 𝑟0 𝑁 (𝑡) ln ( )
𝑑𝑡 𝑁 (𝑡)

Ce modèle repose sur l’hypothèse sous-jacente que la force de mortalité aug-

mente de façon exponentielle avec l’âge (Pletcher, 1999).
A vous de jouer !
• Résolution exacte
• Simulation
• Points d’équilibre et stabilité
• Point d’inflexion
Qu’est-ce qui distingue le modèle de Gompertz du modèle logistique ? Argu-
mentez votre réponse.
Eléments de correction pour le modèle de Gompertz.

2.6.2 Le modèle de von Bertalanffy

En 1934, Karl Ludwig von Bertalanffy a proposé un nouveau modèle de crois-
sance, dont il existe aujourd’hui plusieurs variantes ainsi qu’une version général-
isée livrée en 1959 par F.J. Richards.
Le modèle de von Bertalanffy s’écrit de la manière suivante :

𝑑𝑁 (𝑡)
= 𝑟0 (𝐾 − 𝑁 (𝑡))
𝑑𝑡

A vous de jouer !
• Résolution exacte
• Simulation
• Points d’équilibre et stabilité
• Point d’inflexion
Pour le modèle de von Bertalanffy, posez-vous la question de sa pertinence pour
décrire des données de croissance bactérienne. Argumentez votre réponse.
Eléments de correction pour le modèle de von Bertalanffy.
2.6. VERS D’AUTRES MODÈLES 31

2.6.3 Le modèle de Baranyi

Un modèle très utilisé en microbiologie prédictive parce qu’il permet de décrire

finement la phase de latence de la croissance bactérienne est celui de József
Baranyi, qui publia, avec son collaborateur Terry A. Roberts en 1993, un modèle
de la forme :

𝑑𝑁 (𝑡)
= 𝜇𝛼(𝑡)𝑓(𝑁 (𝑡))𝑁 (𝑡)
𝑑𝑡

dont il donne la solution en 𝑦(𝑡) = ln(𝑁 (𝑡)) (voir ci-dessous). Si complexe soit-
il, ce modèle décrit très bien les différentes phases de la croissance microbienne
– toujours avec une courbe sigmoïde, y compris la phase de latence grâce à la
fonction 𝛼(𝑡), dont la courbe représentative est elle-même en forme de « S ».
xmax <- 10^9
x0 <- 10^3
mumax <- 1.5
lag <- 6
h0 <- mumax*lag
A <- function(x, mumax, h0){
return(x + log(exp(-mumax*x) + exp(-h0) - exp(-mumax*x-h0))/mumax)
}
Baranyi <- function(x, xmax, x0, mumax, lag){
ymax <- log(xmax)
y0 <- log(x0)
return(y0 + mumax*A(x, mumax, lag) - log(1+(exp(mumax*A(x, mumax, lag))-1)/exp(ymax-y0)))
}
par(mar = c(4, 4, 0.2, 0.2))
curve(Baranyi(x, xmax, x0, mumax, lag)/log(2),
from = 0, to = 22, las = 1, lwd = 2,
xlab = "", ylab = "", xaxt = "n", yaxt = "n",
xlim = c(-1, 22))
mtext("Temps", side = 1, line = 0.75)
mtext(expression(Log[2](densité)), side = 2, line = 0.5)
segments(-1, log2(x0), 0, log2(x0), lwd = 2)
 32CHAPTER 2. MODÉLISATION DE LA CROISSANCE DES MICRO-ORGANISMES

Log2(densité)

Temps
L’ ’ ’équation différentielle correspondant au modèle de Baranyi et Roberts est
la suivante :
𝑑𝑁 (𝑡)
= 𝜇𝛼(𝑡)𝑓(𝑁 (𝑡))𝑁 (𝑡)
𝑑𝑡
La solution exprimée en 𝑦(𝑡) = ln(𝑁 (𝑡)) s’écrit en fait :
𝑒𝜇𝐴(𝑡) − 1
𝑦(𝑡) = 𝑦0 + 𝜇𝐴(𝑡) − ln (1 + )
𝑒𝑦𝑚 𝑎𝑥−𝑦0
où la fonction 𝐴(𝑡) s’écrit elle-même comme suit :
1
𝐴(𝑡) = 𝑡 + ln (𝑒−𝜇𝑡 (1 − 𝑒−ℎ0 ) + 𝑒−ℎ0 )
𝜇
Enfin, on a ℎ0 = 𝜇×𝑙𝑎𝑔, où 𝑙𝑎𝑔 est le quatrième paramètre introduit par Baranyi
et Roberts pour inclure la phase de latence dans la modélisation de la croissance
microbienne.

2.7 Comparaison des modèles de croissance bac-

térienne
2.7.1 Comparaison en termes d’accroissement
Pour les trois modèles dont il a été question précédemment (Verhulst, Gompertz,
von Bertalanffy), représentez sur un même graphe les courbes 𝑑𝑁(𝑡)
𝑑𝑡 en fonction
de la variable 𝑁 (𝑡).
Que remarquez-vous ?
2.7. COMPARAISON DES MODÈLES DE CROISSANCE BACTÉRIENNE33

2.7.2 Comparaison des courbes représentatives

Complétez la légende sur le graphe ci-dessous et discutez.

N0 <- 100
K <- 9000
r0 <- 0.85
par(mar=c(4, 4, 0.2, 0.2))
curve(K/(1+(K/N0-1)*exp(-r0*x)), from=0, to=15, las=1, col="black",
lwd=2, xlab="Temps (h)",
ylab="Taille de la population bactérienne (UFC)")
# x représente le temps
abline(h=K/2, lwd=2, lty=2)
curve(K*exp(log(N0/K)*exp(-r0*x)), from=0, to=15, lwd=2,
col="blue", add=TRUE)
abline(h=K/exp(1), lwd=2, lty=2, col="blue")
curve(K-(K-N0)*exp(-r0*x), from=0, to=15, lwd=2,
col="green", add=TRUE)
legend(10,2000, legend=c("....................",
"....................",
"...................."),
lty=1, col=c("black","blue","green"), bty="n",
[Link]="gray", lwd=2)
Taille de la population bactérienne (UFC)

8000

6000

4000

2000
....................
....................
....................
0

0 5 10 15

Temps (h)

Eléments de correction pour la comparaison des modèles.

34CHAPTER 2. MODÉLISATION DE LA CROISSANCE DES MICRO-ORGANISMES

2.8 Calcul du temps de doublement

Sous l’hypothèse du modèle exponentiel, le taux de croissance est constant (𝑟).
On peut alors calculer un temps de doublement de la population bactérienne,
c’est-à-dire un temps 𝜏𝑑 tel que 𝑁 (𝑡 + 𝜏𝑑 ) = 2𝑁 (𝑡). Sachant que 𝑁 (𝑡) = 𝑁0 𝑒𝑟𝑡 ,
il vient :
𝑁 (𝑡 + 𝜏𝑑 ) = 2𝑁 (𝑡)
⇔ 𝑁0 𝑒𝑟(𝑡+𝜏𝑑 ) = 2𝑁0 𝑒𝑟𝑡
⇔ 𝑒𝑟(𝑡+𝜏𝑑 ) = 2𝑒𝑟𝑡
⇔ 𝑒𝑟𝑡 𝑒𝑟𝜏𝑑 = 2𝑒𝑟𝑡
⇔ 𝑒𝑟𝜏𝑑 = 2
⇔ 𝑟𝜏𝑑 = ln 2
⇔ 𝜏𝑑 = ln𝑟2
Le temps de doublement est donc inversement proportionnel au taux de crois-
sance.
Chapter 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 (vari-
able 𝑁 (𝑡)) et de la consommation en substrat (variable 𝑆(𝑡)) à 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 𝑆0 avec un
flux d’entrée 𝐹 , égal au taux de sortie, de sorte que le volume total 𝑉 reste
constant. D’après Edelstein-Keshet (1984).

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

35
36 CHAPTER 3. MODÉLISATION DU CHÉMOSTAT

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 𝐹 , 𝑆0 et 𝑉 .

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 𝑡 dans la 𝑆(𝑡) masse/volume
chambre de culture
Biomasse bactérienne à l’instant 𝑡 dans la chambre 𝑁 (𝑡) CFU/volume
de culture
Concentration en nutriments dans le réservoir 𝑆0 masse/volume
d’entrée
Volume de la chambre de culture 𝑉 volume
Flux d’entrée/sortie 𝐹 volume/temps
Taux de croissance bactérien 𝑟 1/temps
Constante de rendement (cf. § 3.2) 𝑌 = 1/𝛼 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 ordi-
naires 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 :

𝑑𝑁 (𝑡) 𝑁 (𝑡)
= 𝑟𝑁 (𝑡) − 𝐹
𝑑𝑡 𝑉
La partie 𝑟𝑁 (𝑡) correspond à la croissance exponentielle par unité de temps avec
𝑟 (dont l’unité est 1/temps) le taux de croissance malthusien (§ 2.2).
La partie −𝐹 𝑁(𝑡)
𝑉 correspond à ce qui sort de la cuve par unité de temps.
3.3. LA NOTION DE SUBSTRAT LIMITANT 37

La division par 𝑉 assure le respect des dimensions dans l’équation. À cause de

l’unité de 𝐹 , la quantité 𝑉𝐹 s’exprime en [𝑡𝑒𝑚𝑝𝑠]−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 :

𝑟 = 𝑟(𝑆)

Nous verrons ultérieurement quelle hypothèse choisir pour la fonction 𝑟(𝑆) : 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 𝛼 unités de nutriment sont
consommées pour produire 1 unité de bactéries. La quantité 𝑌 = 1/𝛼
représente ainsi le rendement, et on peut écrire :

𝑑𝑆(𝑡) 𝑆(𝑡) 𝑆
= −𝛼𝑟(𝑆)𝑁 (𝑡) − 𝐹 +𝐹 0
𝑑𝑡 𝑉 𝑉
𝑆0
avec −𝛼𝑟(𝑆)𝑁 (𝑡) ce qui est consommé, −𝐹 𝑆(𝑡)𝑉 ce qui sort de la cuve et +𝐹 𝑉
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 lim-
ité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 𝑟(𝑆)) en fonction
de la concentration en nutriments (𝑆) n’est pas proportionnel à 𝑆 mais suit une
relation monotone croissante :

𝑟𝑚𝑎𝑥 𝑆
𝑟(𝑆) =
𝐾𝑠 + 𝑆
38 CHAPTER 3. MODÉLISATION DU CHÉMOSTAT

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 𝑣 à la concentration en substrat [𝑆] ; il inclut des
paramètres constants caractéristiques de l’enzyme :

𝑣𝑚𝑎𝑥 [𝑆]
𝑣([𝑆]) =
𝐾𝑀 + [𝑆]

𝐾𝑀 est la constante de Michaelis spécifique de l’enzyme ; elle a la dimension

d’une concentration (même unité que [𝑆]). Elle est définie par 𝐾𝑀 = 𝑘3𝑘+𝑘2 où
1
les constantes 𝑘𝑖 sont les constantes de vitesse de la réaction chimique suivante
:

𝑘1 𝑘3
𝐸 + 𝑆 ⇌ 𝐸𝑆 → 𝐸 + 𝑃
𝑘2

𝑣 = 𝑘3 × [𝐸𝑆]

𝐾𝑀 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.

0.8

0.6
r(S)

0.4

0.2

0.0

0 2 4 6 8 10

S
3.4. REPARAMÉTRISATION DU MODÈLE 39

Fig. 5 : Représentation graphique du taux de croissance bactérien (𝑟) en fonction

de la concentration en nutriments (𝑆).
• A quoi correspondent mathématiquement les paramètres 𝑟𝑚𝑎𝑥 et 𝐾𝑠 ?
• Quelle est la signification biologique du paramètre 𝑟𝑚𝑎𝑥 ?

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 :
𝑑𝑁(𝑡)
𝑑𝑡 = ( 𝑟𝐾max 𝑆(𝑡)
+𝑆(𝑡) ) 𝑁 (𝑡) −
𝐹
𝑉 𝑁 (𝑡)
{ 𝑠
𝑑𝑆(𝑡)
𝑑𝑡 = −𝛼 ( 𝑟𝐾max 𝑆(𝑡)
)𝑁 (𝑡) − 𝐹
𝑆 (𝑡) + 𝐹
𝑆0
𝑠 +𝑆(𝑡) 𝑉 𝑉

Ce système contient 6 paramètres : 𝑟𝑚𝑎𝑥 , 𝐾𝑠 , 𝐹 , 𝑉 , 𝛼 et 𝑆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 𝑁 (𝑡) = 𝑁 ∗ (𝑡)× 𝑁̂ , 𝑆(𝑡) = 𝑆 ∗ (𝑡)× 𝑆 ̂ et 𝑡 = 𝑡∗ ×𝜏 ,
avec 𝑁 ∗ (𝑡), 𝑆 ∗ (𝑡) et 𝑡∗ de nouvelles variables sans dimension, et 𝑁̂ , 𝑆 ̂ et 𝜏 des
constantes indépendantes du temps. Ainsi :

𝑑(𝑁 ∗ 𝑁)̂ 𝑟max 𝑆 ∗ 𝑆 ̂

⎧ ∗ ̂ 𝐹 ∗ ̂
{ 𝑑(𝑡∗ 𝜏) = ( 𝐾𝑠 +𝑆 ∗ 𝑆 ̂ ) 𝑁 𝑁 − 𝑉 𝑁 𝑁
⎨ ∗ ̂
𝑑(𝑆 𝑆) 𝑟max 𝑆 ∗ 𝑆 ̂ ∗ ̂ 𝐹 ∗ ̂ 𝐹
{ 𝑑(𝑡∗ 𝜏) = −𝛼 ( 𝐾𝑠 +𝑆 ∗ 𝑆 ̂ ) 𝑁 𝑁 − 𝑉 𝑆 𝑆 + 𝑉 𝑆0
⎩

Pour comprendre ce que sont ces nouvelles variables, prenons l’exemple d’une
biomasse bactérienne qui vaut 𝑁 = 105 UFC/ml, alors :

𝑁 = 105 UFC/ml
= 1 (× 105 UFC) /ml
= 100 UFC/𝜇l
= 𝑁 ∗ × 𝑁̂

avec 𝑁 ∗ un nombre sans dimension et 𝑁̂ 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
𝜏 et on divise par 𝑁̂ dans la première équation (resp. par 𝑆 ̂ dans la deuxième
équation), on obtient :

⎧ 𝑑(𝑁 ∗ ) 𝑆∗
{ 𝑑(𝑡∗ ) = 𝜏 𝑟max ( 𝐾 ̂ ∗
) 𝑁 ∗ − 𝜏 𝑉𝐹 𝑁 ∗
𝑠 /𝑆+𝑆
⎨ 𝑑(𝑆 ∗ )
= −𝛼𝜏 𝑟max 𝑁 (
̂ 𝑆∗ ∗
− 𝜏 𝑉𝐹 𝑆 ∗ + 𝜏 𝑉𝐹 𝑆0
{
⎩ 𝑑(𝑡∗ ) 𝑆̂ 𝐾 ̂ ∗)𝑁 𝑆̂
𝑠 𝑆+𝑆
/
40 CHAPTER 3. MODÉLISATION DU CHÉMOSTAT

Si on pose maintenant 𝜏 = 𝑉
𝐹 , 𝑆 ̂ = 𝐾𝑠 et 𝑁̂ = 𝐾𝑠
𝛼𝜏𝑟𝑚𝑎𝑥 , alors il vient :

𝑑𝑁 ∗ 𝑆 ∗
∗ ∗
𝑑𝑡∗ = 𝜏 𝑟max ( 1+𝑆 ∗)𝑁 − 𝑁
{ 𝑑𝑆 ∗ 𝑆∗ ∗ ∗ 𝑆0
𝑑𝑡∗ = − ( 1+𝑆∗ ) 𝑁 − 𝑆 + 𝐾
𝑠

𝑉
• Quelle interprétation peut-on donner au paramètre 𝜏 = 𝐹 ?
𝑆0
Enfin, si on pose 𝛼1 = 𝜏 𝑟𝑚𝑎𝑥 = 𝑉 𝑟𝐹𝑚𝑎𝑥 et 𝛼2 = 𝐾 , alors on obtient un nouveau
𝑠
modèle qui ne contient plus que deux paramètres :
𝑑𝑁 ∗ 𝑆 ∗
∗ ∗
𝑑𝑡∗ = 𝛼1 ( 1+𝑆 ∗)𝑁 − 𝑁
{ 𝑑𝑆 ∗ 𝑆 ∗
∗ ∗
𝑑𝑡∗ = − ( 1+𝑆 ∗ ) 𝑁 − 𝑆 + 𝛼2

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 𝛼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 𝑁 (𝑡) et 𝑆(𝑡) restent constantes toutes les deux simul-
tanément. Mathématiquement, cela revient donc à vérifier :

𝑑𝑁
=0 𝛼 ( 𝑆 )𝑁 − 𝑁 = 0
{ 𝑑𝑡
𝑑𝑆 ⇔ { 1 1+𝑆
𝑆
𝑑𝑡 =0 − ( 1+𝑆 ) 𝑁 − 𝑆 + 𝛼2 = 0

𝑑𝑁 𝑑𝑆
On qualifie de points d’équilibre les solutions du système 𝑑𝑡 = 0 et 𝑑𝑡 = 0.
𝑆2 1
La première équation conduit aux deux solutions 𝑁1 = 0 ou 1+𝑆2 = 𝛼1 , c’est-à-
dire 𝑁1 = 0 et 𝑆2 = 𝛼 1−1 .
1

On considérera dans la suite que 𝛼1 ≠ 1 pour que 𝑆2 existe mathématiquement.

• Si 𝑁1 = 0, alors la deuxième équation du système conduit à 𝑆1 = 𝛼2 . Dans
le plan (𝑁 , 𝑆), on a donc un premier point d’équilibre (𝑁1 , 𝑆1 ) = (0, 𝛼2 ).
– Quelle interprétation biologique pouvez-vous donner à ce point
d’équilibre ?
3.5. ANALYSE MATHÉMATIQUE DU MODÈLE DU CHÉMOSTAT 41

1
• Si 𝑆2 = 𝛼1 −1 , alors la deuxième équation du système conduit à
𝑆2
( 1+𝑆 ) 𝑁2 = 𝛼2 − 𝑆2 , c’est-à-dire 𝑁2 = ( 1+𝑆
𝑆2 ) (𝛼2 − 𝑆2 ) =
2
2
𝛼1 (𝛼2 − 𝑆2 ). On obtient donc un deuxième point d’équilibre
(𝑁2 , 𝑆2 ) = (𝛼1 (𝛼2 − 𝛼 1−1 ), 𝛼 1−1 ).
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 (𝑥𝑒𝑞 , 𝑦𝑒𝑞 ) du
système d’équations :
𝑑𝑥
𝑑𝑡 = 𝑓(𝑥, 𝑦)
{ 𝑑𝑦
𝑑𝑡 = 𝑔(𝑥, 𝑦)

est asymptotiquement stable si et seulement si les deux conditions suivantes

sont réunies :

𝜕𝑓(𝑥𝑒𝑞 ,𝑦𝑒𝑞 ) 𝜕𝑔(𝑥𝑒𝑞 ,𝑦𝑒𝑞 )

𝜕𝑥 + 𝜕𝑦 < 0 (critère − 1)
{ 𝜕𝑓(𝑥𝑒𝑞 ,𝑦𝑒𝑞 ) 𝜕𝑔(𝑥𝑒𝑞 ,𝑦𝑒𝑞 ) 𝜕𝑓(𝑥𝑒𝑞 ,𝑦𝑒𝑞 ) 𝜕𝑔(𝑥𝑒𝑞 ,𝑦𝑒𝑞 )
𝜕𝑥 . 𝜕𝑦 − 𝜕𝑦 . 𝜕𝑥 >0 (critère − 2)

𝜕𝑓(𝑥 ,𝑦 ) 𝜕𝑓(𝑥 ,𝑦 ) 𝜕𝑔(𝑥 ,𝑦 ) 𝜕𝑔(𝑥 ,𝑦 )

où les termes 𝑒𝑞 𝑒𝑞
𝜕𝑥 , 𝑒𝑞 𝑒𝑞
𝜕𝑦 , 𝑒𝑞 𝑒𝑞
𝜕𝑥 , et 𝑒𝑞 𝑒𝑞
𝜕𝑦 sont les dérivées
partielles des fonctions 𝑓(𝑥, 𝑦) et 𝑔(𝑥, 𝑦) par rapport aux variables 𝑥 et 𝑦, re-
spectivement, évaluées en 𝑥 = 𝑥𝑒𝑞 et 𝑦 = 𝑦𝑒𝑞 .
Dans le cas du modèle du chémostat, nous avons :
𝑆
𝑓 (𝑁 , 𝑆) = 𝛼1 ( 1+𝑆 )𝑁 − 𝑁
𝑆
𝑔 (𝑁 , 𝑆) = − ( 1+𝑆 ) 𝑁 − 𝑆 + 𝛼2

𝑥 ′ 1
• Vérifiez que ( 1+𝑥 ) = (1+𝑥)2

Pour le modèle du chémostat, le calcul des dérivées partielles conduit à :

𝜕𝑓(𝑁,𝑆) 𝑆
𝜕𝑁 = 𝛼1 ( 1+𝑆 )−1
𝜕𝑓(𝑁,𝑆) 𝑁
𝜕𝑆 = 𝛼1 2
(1+𝑆)
𝜕𝑔(𝑁,𝑆) 𝑆
𝜕𝑁 = − ( 1+𝑆 )
𝜕𝑔(𝑁,𝑆) 𝑁
𝜕𝑆 =− 2 − 1
(1+𝑆)
42 CHAPTER 3. MODÉLISATION DU CHÉMOSTAT

Il nous faut maintenant calculer ces dérivés partielles en chacun des points
d’équilibre (𝑁1 , 𝑆1 ) et (𝑁2 , 𝑆2 ). Pour simplifier les calculs, posons :
𝛼2
𝐴=
1 + 𝛼2
et
𝑁2
𝐵=
(1 + 𝑆2 )2
Les deux quantités 𝐴 et 𝐵 sont > 0.
• Au point d’équilibre (𝑁1 , 𝑆1 ) = (0, 𝛼2 ), il vient :
𝜕𝑓(𝑁1 ,𝑆1 ) 𝛼2
𝜕𝑁 = 𝛼1 ( 1+𝛼 ) − 1 = 𝛼1 𝐴 − 1
2
𝜕𝑓(𝑁1 ,𝑆1 )
𝜕𝑆 =0
𝜕𝑔(𝑁1 ,𝑆1 ) 𝛼2
𝜕𝑁 = − ( 1+𝛼 ) = −𝐴
2
𝜕𝑔(𝑁1 ,𝑆1 )
𝜕𝑆 = −1

Ainsi :
𝜕𝑓(𝑁1 ,𝑆1 )
𝜕𝑁 + 𝜕𝑔(𝑁𝜕𝑆1 ,𝑆1 ) = 𝛼1 𝐴 − 2
𝜕𝑓(𝑁1 ,𝑆1 ) 𝜕𝑔(𝑁1 ,𝑆1 )
𝜕𝑁 . 𝜕𝑆 − 𝜕𝑓(𝑁 1 ,𝑆1 ) 𝜕𝑔(𝑁1 ,𝑆1 )
𝜕𝑆 . 𝜕𝑁 = 1 − 𝛼1 𝐴

Comme nous l’avons vu précédemment, le deuxième point d’équilibre a du sens

biologiquement si 𝛼1 > 1 et 𝛼2 (𝛼1 − 1) > 1. Nous voyons donc que si ces deux
conditions sont vérifiées, alors :

𝛼2 1 + 𝛼2 − 𝛼1 𝛼2 1 − 𝛼2 (𝛼1 − 1)
1 − 𝛼 1 𝐴 = 1 − 𝛼1 ( )= = <0
1 + 𝛼2 1 + 𝛼2 1 + 𝛼2

Par conséquent, le point d’équilibre (𝑁1 , 𝑆1 ) est instable.

• Le deuxième point d’équilibre (𝑁2 , 𝑆2 ) s’écrit (𝛼1 (𝛼2 − 𝛼 1−1 ), 𝛼 1−1 ). Il n’a
1 1
𝑆2 1
de sens biologique que si 𝛼1 > 1 et 𝛼2 (𝛼1 − 1) > 1, et il vérifie 1+𝑆2 = 𝛼1 .
Ainsi, il vient :
𝜕𝑓(𝑁2 ,𝑆2 )
𝜕𝑁 =0
𝜕𝑓(𝑁2 ,𝑆2 )
𝜕𝑆 = 𝛼1 𝐵
𝜕𝑔(𝑁2 ,𝑆2 )
𝜕𝑁 = − 𝛼1
1
𝜕𝑔(𝑁2 ,𝑆2 )
𝜕𝑆 = −𝐵 − 1
Le calcul des conditions de stabilité conduit à :
𝜕𝑓(𝑁2 ,𝑆2 )
𝜕𝑁 + 𝜕𝑔(𝑁𝜕𝑆2 ,𝑆2 ) = −𝐵 − 1 < 0
𝜕𝑓(𝑁2 ,𝑆2 ) 𝜕𝑔(𝑁2 ,𝑆2 )
𝜕𝑁 . 𝜕𝑆 − 𝜕𝑓(𝑁 2 ,𝑆2 ) 𝜕𝑔(𝑁2 ,𝑆2 )
𝜕𝑆 . 𝜕𝑁 =𝐵>0

Par conséquent, le point d’équilibre (𝑁2 , 𝑆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.6. REPRÉSENTATIONS GRAPHIQUES DES SOLUTIONS D’UN SYSTÈME DYNAMIQUE43

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 (𝑁2 , 𝑆2 ) qui fait sens biologiquement
si 𝛼1 > 1 et 𝛼2 (𝛼1 − 1) > 1 et qui est asymptotiquement stable.

• La condition 𝛼1 > 1 est équivalente à 𝑟𝑚𝑎𝑥 > 𝜏1 ⇔ 𝑟 1 < 𝑉𝐹 . 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 𝜏𝑑 = 𝑟ln 2 .
𝑚𝑎𝑥
Ainsi, la condition 𝛼1 > 1 équivaut à 𝜏𝑑 < ln 2 𝑉𝐹 , c’est-à-dire qu’elle
impose que 𝜏𝑑 soit inférieur au temps nécessaire pour vider la chambre
de culture (× 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.

𝑆0 1
• La seconde condition 𝛼2 (𝛼1 − 1) > 1 équivaut à 𝐾𝑠 > 𝛼1 −1 ,
c’est-à-dire
𝑆0
> 𝑆2∗ , soit encore 𝑆0 > 𝐾𝑠 𝑆2∗ . Si on ̂
se rappelle que 𝑆 = 𝐾𝑠 , il vient
𝐾𝑠
̂ ∗ . La quantité 𝑆𝑆
𝑆0 > 𝑆𝑆 ̂ ∗ correspondà la valeur de 𝑆2 exprimée dans
2 2
l’unité de départ avant le re-dimensionnement, c’est-à-dire exprimée en
masse/volume, comme l’est 𝑆0 . En conclusion, nous tombons sur un ré-
sultat trivial : la concentration en nutriments à l’équilibre (𝑆𝑆 ̂ ∗)
2
ne peut pas être supérieure à la concentration du stock de nu-
triment (𝑆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

𝑑𝑁
𝑑𝑡 = 𝑓(𝑁 , 𝑆)
{ 𝑑𝑆
𝑑𝑡 = 𝑔(𝑁 , 𝑆)

les solutions sont de la forme (𝑁 (𝑡), 𝑆(𝑡)).

Ces solutions peuvent être représentées graphiquement soit dans un plan (𝑁 , 𝑆),
qu’on appelle le plan de phase, soit dans des plans (𝑡, 𝑁 (𝑡)) ou (𝑡, 𝑆(𝑡)). Dans
le plan de phase, la représentation graphique d’une solution (𝑁 (𝑡), 𝑆(𝑡)) des-
sine une trajectoire depuis une condition initiale (𝑁 (0), 𝑆(0)) et par exemple
jusqu’à un point d’équilibre asymptotiquement stable (𝑁𝑒𝑞 , 𝑆𝑒𝑞 ) (Fig. 6).
 44 CHAPTER 3. MODÉLISATION DU CHÉMOSTAT

2.5

(N(0),S(0))
2.0

1.5
v
S(t)

1.0
(Neq,Seq)

0.5

0.0

0.0 0.5 1.0 1.5 2.0 2.5

N(t)

Fig. 6 : Représentation graphique d’une solution (𝑁 (𝑡), 𝑆(𝑡)) dans le plan de

phase depuis une condition initiale (𝑁 (0), 𝑆(0)) jusqu’à un point d’équilibre
asymptotiquement stable (𝑁𝑒𝑞 , 𝑆𝑒𝑞 ).

Sans résoudre analytiquement le système dynamique, la solution (𝑁 (𝑡), 𝑆(𝑡))

est en fait connue indirectement par les fonctions 𝑓(𝑁 , 𝑆) et 𝑔(𝑁 , 𝑆) qui carac-
térisent en chaque point (et donc à chaque pas de temps) le vecteur vitesse
:

𝑑𝑁
V=( 𝑑𝑡 )
𝑑𝑆
𝑑𝑡

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).
3.6. REPRÉSENTATIONS GRAPHIQUES DES SOLUTIONS D’UN SYSTÈME DYNAMIQUE45

2.5

2.0

1.5
S(t)

1.0

0.5

0.0

0.0 0.5 1.0 1.5 2.0 2.5

N(t)

𝑑𝑁 𝑑𝑆
Fig. 7 : Champ de vecteurs du système défini par 𝑑𝑡 = 𝑓(𝑁 , 𝑆) et 𝑑𝑡 = 𝑔(𝑁 , 𝑆).

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

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 :
𝑑𝑁
V=( 𝑑𝑡 )
0

Par conséquent, les isoclines nulles horizontales sont les courbes solution
de l’équation :
𝑑𝑆
= 0 ⇔ 𝑔(𝑁 , 𝑆) = 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 :

0
V=( 𝑑𝑆 )
𝑑𝑡
46 CHAPTER 3. MODÉLISATION DU CHÉMOSTAT

Par conséquent, les isoclines nulles horizontales sont les courbes solution
de l’équation :
𝑑𝑁
= 0 ⇔ 𝑓(𝑁 , 𝑆) = 0
𝑑𝑡
Ainsi, on peut compléter le champ de vecteurs avec les isoclines nulles dans le
plan de phase (Fig. 8).

2.5

2.0

1.5
S(t)

1.0

0.5

0.0

0.0 0.5 1.0 1.5 2.0 2.5

N(t)

Fig. 8 : Représentation des isoclines nulles horizontales (en bleu) et verticales

(en vert) du système défini par 𝑑𝑁 𝑑𝑆
𝑑𝑡 = 𝑓(𝑁 , 𝑆) et 𝑑𝑡 = 𝑔(𝑁 , 𝑆).

Puisque les points d’équilibre vérifient simultanément 𝑑𝑁 𝑑𝑆

𝑑𝑡 = 0 et 𝑑𝑡 = 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)
3.6. REPRÉSENTATIONS GRAPHIQUES DES SOLUTIONS D’UN SYSTÈME DYNAMIQUE47

}
# 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,
[Link] = 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)
48 CHAPTER 3. MODÉLISATION DU CHÉMOSTAT

2.5

2.0

1.5
S(t)

1.0

0.5

0.0

0.0 0.5 1.0 1.5 2.0 2.5

N(t)

𝑑𝑁 𝑑𝑆
Fig. 9 : Trajectoires du système du système défini par 𝑑𝑡 = 𝑓(𝑁 , 𝑆) et 𝑑𝑡 =
𝑔(𝑁 , 𝑆).

3.6.3 Les chroniques

Comme évoqué précédemment, on peut aussi représenter les solutions d’un sys-
tème dynamique dans les plans (𝑡, 𝑁 (𝑡)) et (𝑡, 𝑆(𝑡)). On appelle ces représenta-
tions 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), [Link]=FALSE,
xlab=expression(t),ylab="")
legend("topright", legend=c(expression(N(t)), expression(S(t))),
col=c("orange","gray"), lwd=2, bty="n")
3.6. REPRÉSENTATIONS GRAPHIQUES DES SOLUTIONS D’UN SYSTÈME DYNAMIQUE49

2.5 N(t)
S(t)

2.0

1.5

1.0

0.5

0.0

0 5 10 15 20 25 30

t
𝑑𝑁
Fig. 10 : Chroniques associées au système défini par par 𝑑𝑡 = 𝑓(𝑁 , 𝑆) et
𝑑𝑆
𝑑𝑡 = 𝑔(𝑁 , 𝑆).

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="[Link]", 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 ap-
proximation numérique grâce à la fonction findEquilibrium() ; cela nécessite
50 CHAPTER 3. MODÉLISATION DU CHÉMOSTAT

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 :
𝑑𝑁 𝑆
𝑑𝑡 = 𝛼1 ( 1+𝑆 )𝑁 − 𝑁
{ 𝑑𝑆 𝑆
𝑑𝑡 = − ( 1+𝑆 ) 𝑁 − 𝑆 + 𝛼2

Il faut en premier lieu choisir des valeurs de paramètres pour 𝛼1 et 𝛼2 , sans

oublier de respecter les deux conditions : 𝛼1 > 1 et 𝛼2 (𝛼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 (𝛼1 , 𝛼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) :
𝑑𝑁(𝑡) 𝑟max 𝑆(𝑡) 𝐹
𝑑𝑡 = ( 𝐾𝑠 +𝑆(𝑡) ) 𝑁 (𝑡) − 𝑉 𝑁 (𝑡)
{ 𝑑𝑆(𝑡) 𝑟max 𝑆(𝑡) 𝐹 𝐹
𝑑𝑡 = −𝛼 ( 𝐾 +𝑆(𝑡) ) 𝑁 (𝑡) − 𝑉 𝑆 (𝑡) + 𝑉 𝑆0
𝑠
3.7. RETOUR AU MODÈLE DU CHÉMOSTAT 51

• 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 com-
portement 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.
52 CHAPTER 3. MODÉLISATION DU CHÉMOSTAT
Chapter 4

Vitesse d’une réaction

enzymatique

Dans cette partie, nous allons revenir sur le modèle de Michaelis-Menten (§

3.3) proposé en 1913 par l’allemand Leonor Michaelis et la canadienne Maud
Menten pour décrire la vitesse (𝑣) d’une réaction chimique catalysée par une
enzyme agissant sur un substrat ([𝑆]) unique pour donner un produit (Michaelis
et Menten, 1913) :

𝑣𝑚𝑎𝑥 [𝑆]
𝑣([𝑆]) =
𝐾𝑀 + [𝑆]

On parle de cinétique enzymatique. Il s’agit d’une relation monotone crois-

sante sans point d’inflexion, telle que lim 𝑣 ([𝑆]) = 𝑣max et 𝑣(𝐾𝑀 ) = 𝑣𝑚𝑎𝑥
2 .
[𝑆]→+∞

Dans cette relation, [𝑆] est la variable indépendante explicative (on

l’appelle aussi variable de contrôle), tandis que 𝑣 est la variable dépen-
dante à expliquer.

4.1 Les observations

Nous allons travailler avec des données tirées de la littérature scientifique (Fig.
11).

53
54 CHAPTER 4. VITESSE D’UNE RÉACTION ENZYMATIQUE

Jeu de données n°1 Jeu de données n°2 Jeu de données n°3

0.5 14 12

12
10
0.4

10
8
0.3
8
Vitesse

Vitesse

Vitesse
6
6
0.2
4
4

0.1
2
2

0.0 0 0

0.0 0.5 1.0 1.5 0 20 40 60 80 0 50 100 150

Concentration en substrat Concentration en substrat Concentration en substrat

Fig. 11: Représentation graphique de la relation vitesse versus concentration

en susbtrat pour trois jeux de données tirés de la littérature (Ritchie et Prvan,
1996).
La tendance moyenne des données semble suivre le modèle de Michaelis-Menten
avec, pour chaque jeu de données, des valeurs des paramètres 𝑣𝑚𝑎𝑥 et 𝐾𝑀
différentes.
Télécharger le jeu de données n°1
Télécharger le jeu de données n°2
Télécharger le jeu de données n°3

4.2 La méthode de régression

4.2.1 Généralités
Soit (𝑥𝑖 , 𝑦𝑖 ) les couples de valeurs observées pour [𝑆] et 𝑣 dans un jeu de données,
avec 𝑖 = 1, ..., 𝑛 et 𝑛 le nombre d’observations dans le jeu de données.
Si on suppose que la concentration en substrat peut être déterminée avec suff-
isamment de précision, alors on peut considérer que les 𝑥𝑖 sont connus “sans
erreur”. Par contre, il est souvent difficile expérimentalement de mesurer la
vitesse d’une réaction chimique, ce qui nous conduit à considérer que les 𝑦𝑖
seront distribuées aléatoirement et selon une certaine loi autour de la tendance
moyenne décrite par le modèle de Michaelis-Menten. La variable 𝑣 étant une
variable quantitative continue, il est raisonnable de supposer qu’elle se distribue
4.2. LA MÉTHODE DE RÉGRESSION 55

selon une loi de probabilité gaussienne, c’est-à-dire une loi normale (ou loi de
Gauss). Ainsi :
𝑦𝑖 ∼ 𝒩(𝑣(𝑥𝑖 ), 𝜎) ⇔ 𝑦𝑖 = 𝑣(𝑥𝑖 ) + 𝜀𝑖

avec 𝜀𝑖 ∼ 𝒩(0, 𝜎). On appelle les 𝜀𝑖 les résidus, c’est-à-dire les écarts entre
les valeurs observées et les valeurs théoriques moyennes prédites par la partie
déterministe du modèle.
Il s’agit d’un modèle statistique dans lequel la fonction 𝑣(𝑥𝑖 ) correspond à la
partie déterministe et la loi normale à la partie stochastique. C’est un modèle qui
compte trois paramètres, 𝑣𝑚𝑎𝑥 , 𝐾𝑀 et 𝜎, dont nous allons chercher à estimer les
valeurs pour chacun des jeux de données ci-dessus. Ainsi, nous allons chercher à
optimiser la valeur des paramètres pour que la courbe représentant la tendance
moyenne passe “au mieux” entre les points ; nous allons faire ce que l’on appelle
de la régression.
Le terme “régression” a été introduit par Francis Galton à la suite d’une étude
sur la taille des descendants de personnes de grande taille, qui diminue de généra-
tions en générations vers une taille moyenne, c’est-à-dire que leur taille régresse
(Galton, 1890).

4.2.2 Le critère des moindres carrés

Lorsque l’on a affaire à un modèle statistique gaussien, l’optimation des
paramètres consiste à minimiser la Somme des Carrés des Écarts (𝑆𝐶𝐸) entre
les observations et la valeur moyenne prédite par le modèle, c’est-à-dire à
minimiser la fonction suivante :

𝑛
2
𝑆𝐶𝐸 (Θ) = ∑ (𝑦𝑖 − 𝑣 (𝑥𝑖 ))
𝑖=1

avec Θ = (𝑣𝑚𝑎𝑥 , 𝐾𝑀 ) le vecteur des paramètres à estimer dans le cas du modèle

de Michaelis-Menten.
On parle d’estimation selon le critère des moindres carrés.
Pour être exhaustif, il y a un troisième paramètre à estimer, 𝜎, lié à la distri-
bution normale supposée des observations autour de v(𝑥𝑖 ). Celui-ci est estimé
une fois le processus de minimisation terminé par :

𝑆𝐶𝐸min
𝜎̂ = √
𝑛−2

Graphiquement, minimiser la 𝑆𝐶𝐸 revient à minimiser la somme des distances

représentées par les segments noirs sur la figure ci-dessous (Fig. 12).
56 CHAPTER 4. VITESSE D’UNE RÉACTION ENZYMATIQUE

Jeu de données n°1

0.7

0.6

0.5
Vitesse

0.4

0.3

0.2

0.1

0.0

0.0 0.5 1.0 1.5

Concentration en substrat

Fig. 12 : Régression du modèle de Michaelis-Menten sur le jeu de données

numéro 1.
Dans le cas du modèle de Michaelis-Menten :
• la partie déterministe est non linéaire, c’est-à-dire que les dérivées par-
tielles de 𝑣([𝑆]) par rapport aux deux paramètres 𝑣𝑚𝑎𝑥 et 𝐾𝑀 dépendent
encore des paramètres ;
• la partie stochastique est gaussienne.
On parle donc de régression non linéaire gaussienne et c’est ce qui nous
occupera pour cette troisième partie. Néanmoins, il existe d’autres formes de ré-
gression : la régression linéaire gaussienne, la régression linéaire non gaussienne
(par exemple, la régression logistique), la régression non linéaire non gaussi-
enne (par exemple, un modèle concentration-effet pour des données de survie
binaires).

4.2.3 Algorithmes de minimisation

La plupart du temps, en régression non linéaire gaussienne, la minimisation de
la 𝑆𝐶𝐸 n’admet pas de solution analytique, c’est-à-dire pas de solution exacte.
Cette minimisation se fait donc par résolution numérique selon une méthode
itérative qui impose :
• une évaluation initiale “à l’oeil” des paramètres ;
• des calculs itératifs pour converger vers le vecteur des paramètres qui
minimise la SCE.
4.3. TRANSFORMATIONS LINÉAIRES 57

Parmi les algorithmes de minimisation les plus connus, on trouve : l’algorithme

de Gauss-Newton, l’algorithme de Golub-Pereyra ou encore l’algorithme de
Levenberg-Marquardt. Ces algorithmes reposent tous sur des déplacements
à petits pas dans la direction de plus grande pente de la fonction SCE, en
linéarisant le modèle à chaque itération ; ils sont pour la plupart implémentés
dans la fonction nls du logiciel R.

A l’époque de Michaelis-Menten, et pendant de nombreuses années ensuite,

l’utilisation de ces algorithmes n’était pas permise, soit parce que les méthodes
n’avaient pas encore été découvertes, soit parce qu’elles nécessitaient de gros
moyens de calculs. Il a donc fallu ruser ; en particulier, plusieurs méthodes ont
été proposées pour transformer les données de cinétique enzymatique de façon
à produire une relation linéaire et ainsi faciliter l’estimation des paramètres
(Atkins et Nimmo, 1975).

4.3 Transformations linéaires

S’il existe près d’une dizaine de transformations linéaires, nous ne verrons que
les trois transformations les plus couramment utilisées (Ritchie et Prvan, 1996)
. La rédaction de ce paragraphe est très largement inspirée de la fiche tdr47
élaborée par J.R. Lobry.

4.3.1 La transformation de Lineweaver & Burk

La transformation linéaire de Lineweaver & Burk (1934) consiste à représenter

𝑦 = 1/𝑣 en fonction de 𝑥 = 1/[𝑆] (Fig. 13) :

1 𝐾 1 1
= 𝑀 +
𝑣 𝑣𝑚𝑎𝑥 [𝑆] 𝑣𝑚𝑎𝑥

Il s’agit bien d’une équation de droite dont la pente est donnée par 𝑣𝐾𝑀 et
𝑚𝑎𝑥
l’ordonnée à l’origine par 1/𝑣𝑚𝑎𝑥 . A cela s’ajoute le paramètre 𝜎 de la loi
normale.
58 CHAPTER 4. VITESSE D’UNE RÉACTION ENZYMATIQUE

Données brutes Transformation de Lineweaver & Burk

0.5 8

0.4
6

0.3

1/v
4
v

0.2

2
0.1

0.0 0

0.0 0.5 1.0 1.5 0 2 4 6 8

[S] 1/[S]

Fig. 13 : Comparaison des données brutes et transformées selon Lineweaver &

Burk sur le jeu de données n°1.

4.3.2 La transformation de Eadie et Hofstee

La transformation linéaire de Eadie et Hofstee (Eadie 1942, Hofstee 1952) con-

siste à représenter 𝑦 = 𝑣 en fonction de 𝑥 = 𝑣/[𝑆] (Fig. 14) :

𝑣
𝑣([𝑆]) = −𝐾𝑀 + 𝑣𝑚𝑎𝑥
[𝑆]

Il s’agit bien d’une équation de droite dont la pente est donnée par −𝐾𝑀 (< 0)
et l’ordonnée à l’origine par 𝑣𝑚𝑎𝑥 . A cela s’ajoute le paramètre 𝜎 de la loi
normale.
4.3. TRANSFORMATIONS LINÉAIRES 59

Données brutes Transformation de Eadie et Hofstee

0.5 0.5

0.4 0.4

0.3 0.3
v

v
0.2 0.2

0.1 0.1

0.0 0.0

0.0 0.5 1.0 1.5 0.0 0.4 0.8 1.2

[S] v/[S]

Fig. 14 : Comparaison des données brutes et transformées selon Eadie et Hofstee

sur le jeu de données n°1.

4.3.3 La transformation de Wilkinson

La transformation linéaire de Wilkinson (1961) consiste à représenter 𝑦 = [𝑆]/𝑣

en fonction de 𝑥 = [𝑆] (Fig. 15) :

[𝑆] 1 𝐾
= [𝑆] + 𝑀
𝑣 𝑣𝑚𝑎𝑥 𝑣𝑚𝑎𝑥

Il s’agit bien d’une équation de droite dont la pente est donnée par 1/𝑣𝑚𝑎𝑥 et
l’ordonnée à l’origine par 𝐾𝑀 /𝑣𝑚𝑎𝑥 . A cela s’ajoute le paramètre 𝜎 de la loi
normale.
60 CHAPTER 4. VITESSE D’UNE RÉACTION ENZYMATIQUE

Données brutes Transformation de Wilkinson

0.5 3.0

2.5
0.4

2.0
0.3

[S]/v
1.5
v

0.2
1.0

0.1
0.5

0.0 0.0

0.0 0.5 1.0 1.5 0.0 0.5 1.0 1.5

[S] [S]

Fig. 15 : Comparaison des données brutes et transformées selon Wilkinson sur

le jeu de données n°1.

4.4 Retour à la régression gaussienne

Lorsque l’on fait de la régression, on cherche donc à estimer les paramètres d’un
modèle à partir de données expérimentales ; on dit qu’on cherche à ajuster un
modèle à des données. Nous allons maintenant voir comment procéder grâce au
logiciel R (pour plus de détails, voir ici). Nous utiliserons le jeu de données
numéro 1 pour illustrer les différents points.

4.4.1 Ajuster un modèle linéaire gaussien avec R

Nous travaillerons avec la transformation linéaire de Lineweaver & Burk dont
la partie déterministe s’écrit :
𝐾𝑀 1
𝑦= 𝑥+
𝑣𝑚𝑎𝑥 𝑣𝑚𝑎𝑥
avec 𝑥 = 1/[𝑆] et 𝑦 = 1/𝑣.
1. On commence par importer les données sous la forme d’un tableau à
deux colonnes : S pour la concentration substrat (notre variable explicative) et
V pour la vitesse de la réaction enzymatique (notre variable à expliquer).
# Création du tableau de données
dataMM1 <- [Link]("[Link]", header = TRUE)
# ATTENTION au chemin d'accès à votre propre fichier "[Link]"
4.4. RETOUR À LA RÉGRESSION GAUSSIENNE 61

# Visualisez le jeu de données

par(mar = c(4, 4, 0.2, 0.2)) # réglage des marges
plot(dataMM1, las = 1, pch = 19)

0.50

0.45

0.40

0.35
V

0.30

0.25

0.20

0.15

0.2 0.4 0.6 0.8 1.0 1.2 1.4

S
n <- nrow(dataMM1) # Nombres d'observations dans le jeu de données
# Transformez les données selon Lineweaver \& Burk
x <- 1/dataMM1$S
y <- 1/dataMM1$V

2. On ajuste le modèle linéaire de Lineweaver & Burk à l’aide de la fonction

lm (pour “linear model”).
m <- lm(y ~ x)
# On demande les résultats de l'ajustement
summary(m)

Call:
lm(formula = y ~ x)

Residuals:
1 2 3 4 5 6
-0.40673 0.71770 -0.02187 0.03369 -0.12712 -0.19567

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.70846 0.30327 5.634 0.00489 **
62 CHAPTER 4. VITESSE D’UNE RÉACTION ENZYMATIQUE

x 0.75279 0.07815 9.632 0.00065 ***

---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4291 on 4 degrees of freedom

Multiple R-squared: 0.9587, Adjusted R-squared: 0.9483
F-statistic: 92.78 on 1 and 4 DF, p-value: 0.0006497
On obtient en sortie les valeurs des résidus (Residuals), c’est-à-dire les écarts
entre observations et valeurs calculées par le modèle, les estimations moyennes
(Estimate) des deux paramètres 1/𝑣𝑚𝑎𝑥 (Intercept) et 𝐾𝑀 /𝑣𝑚𝑎𝑥 (x), avec leur
écart-type (Std. Error). L’estimation moyenne de 𝜎 est donnée par la quantité
Residual standard error (ici 0.4291).
# Explorer l'objet 'm' à l'aide de la fonction names()
names(m)

[1] "coefficients" "residuals" "effects" "rank"

[5] "[Link]" "assign" "qr" "[Link]"
[9] "xlevels" "call" "terms" "model"
# Vous constatez que vous avez accès à différents sous-objets.
# Par exemple les coefficients de la droite de régression
m$coefficients

(Intercept) x
1.7084585 0.7527938
Précision et erreur relative
On peut se faire une idée de la précision avec laquelle un paramètre 𝛼 a été
estimé en évaluant l’erreur relative comme suit :

𝑆𝐸(𝛼)
𝐸𝑅(𝛼) = 𝑡
𝑚𝑜𝑦(𝛼) (0.975;𝑛−𝑝)

avec 𝑚𝑜𝑦(𝛼) la moyenne du paramètre 𝛼, 𝑆𝐸(𝛼) l’écart-type (ou erreur stan-

dard, standard error en anglais) du paramètre 𝛼, 𝑝 le nombre de paramètres
dans la partie déterministe du modèle (ici 𝑝 = 2), 𝑛 la taille du jeu de données
(ici 𝑛 = 6) et 𝑡(0.975;𝑛−𝑝) le quantile 97.5% de la loi de Student à 𝑛 − 𝑝 degrés
de liberté (𝑡(0.975;𝑛−𝑝) ∼ 2).
La précision (exprimée en %) est jugée d’autant meilleure que 𝐸𝑅 est petite.
estim <- summary(m)$coefficients[,"Estimate"] # estimations moyennes
SE <- summary(m)$coefficients[,"Std. Error"] # écart-types
ER <- 100*qt(0.975, n-2)*SE/estim # erreurs relatives
ER # affichage (en %)

(Intercept) x
4.4. RETOUR À LA RÉGRESSION GAUSSIENNE 63

49.28430 28.82484
• Qu’en pensez-vous ?
3. On représente ensuite les données sur lesquelles on peut superposer
la droite ajustée, fonction abline, et éventuellement aussi les résidus que l’on
récupère dans l’objet fitted(m).
par(mfrow=c(1,1)) # pour n'avoir qu'un seul graphe dans la fenêtre graphique
par(mar=c(4, 4, 0.2, 0.2)) # réglage des marges
plot(x, y, pch = 19)
abline(m) # superposition de la droite de régression
segments(x, y, x, fitted(m)) # ajout des distances entre observations et valeurs théoriques
6
5
y

4
3
2

1 2 3 4 5 6 7

• Améliorez l’esthétique du graphe si le coeur vous en dit.

On fait ensuite l’analyse des résidus.
En effet, notre modèle statistique complet s’écrit :

𝑌𝑜𝑏𝑠 = 𝑦𝑡ℎ𝑒𝑜 + 𝜀

avec 𝜀 ∼ 𝒩(0, 𝜎). Autrement dit, les résidus doivent vérifier plusieurs hy-
pothèses :
• ils doivent être indépendants les uns des autres : on parle d’absence
d’auto-corrélation entre les résidus ;
64 CHAPTER 4. VITESSE D’UNE RÉACTION ENZYMATIQUE

• ils doivent suivre une loi normale ;

• la variance de la loi normale des résidus est constante et égale à 𝜎2 : on

parle d’homoscédasticité.

4. On représente graphiquement les résidus

# Graphe des résidus en fonction des valeurs prédites
sigma <- summary(m)$sigma # écart-type résiduel
# fitted(m) fournit les valeurs théoriques
# residuals(m) fournit les résidus
plot(fitted(m), residuals(m), pch = 19,
main="Résidus en fonction des valeurs prédites",
las=1, ylim=c(-2.2*sigma,2.2*sigma))
abline(h=0, col="blue") # droite horizontale d'abcisse 0
# On trace ensuite les droites horizontales pour -2*sigma et +2*sigma
abline(h=c(-2*sigma,2*sigma), lty=2, col="blue") # pointillés

Résidus en fonction des valeurs prédites

1.0

0.5
residuals(m)

0.0

−0.5

−1.0
3 4 5 6 7

fitted(m)

Sur ce type de graphe, on s’attend à une répartition aléatoire des résidus

selon une loi normale de variance 𝜎2 constante avec environ 95% des résidus dans
l’intervalle [−2𝜎; 2𝜎], c’est-à-dire entre les deux droites horizontales en pointillés
bleus.

Que faut-il repérer sur le graphe des résidus ?

4.4. RETOUR À LA RÉGRESSION GAUSSIENNE 65

Les résidus ne sont pas distribués aléatoirement → on rejette le modèle

On voit un graphe en entonnoir (hétéroscédasticité) → on rejette le modèle

Il y a des résidus extrêmes (outliers) → on rejette le modèle

66 CHAPTER 4. VITESSE D’UNE RÉACTION ENZYMATIQUE

→ Rien ne permet de rejeter le modèle.

5. On vérifie la normalité des résidus

On peut représenter graphiquement leurs quantiles en fonction des quantiles

d’une loi normale. On s’attend sur ce graphe à ce que les points s’alignent le
long de droite de Henry (en trait plein ci-dessous).
# Graphe des quantiles de nos résidus (axes des x)
# en fonction des quantiles théoriques (axes des y)
# sous l'hypothèse de la loi normale N(0, sigma)
qqnorm(residuals(m), las = 1, pch = 19)
# Ajout de la droite de Henry
qqline(residuals(m))

Normal Q−Q Plot

0.6
Sample Quantiles

0.4

0.2

0.0

−0.2

−0.4

−1.0 −0.5 0.0 0.5 1.0

Theoretical Quantiles

Moyennes et variances
4.4. RETOUR À LA RÉGRESSION GAUSSIENNE 67

Il faut veiller au fait qu’on travaille ici avec la transformation de Lineweaver

& Burk dont les paramètres sont 1/𝑣𝑚𝑎𝑥 et 𝐾𝑀 /𝑣𝑚𝑎𝑥 . Ce sont donc ces deux
paramètres que l’on estime et non pas directement les paramètres 𝑣𝑚𝑎𝑥 et 𝐾𝑀
du modèle de Michaelis-Menten.
De manière générale, on retiendra les résultats suivants :
• Relation entre les moyennes :
1 1
𝑚𝑜𝑦 ( ) =
𝛼 𝑚𝑜𝑦(𝛼)

𝛼 𝑚𝑜𝑦(𝛼)
𝑚𝑜𝑦 ( ) =
𝛽 𝑚𝑜𝑦(𝛽)

• Une approximation de l’erreur standard (𝑆𝐸) sur l’inverse de 𝛼 est :

1 𝑆𝐸(𝛼)
𝑆𝐸 ( ) ∼
𝛼 𝛼2

• Si on néglige les termes de covariance, une approximation de l’erreur stan-

dard d’un rapport est :

2 2
𝛼 𝛼 𝑆𝐸(𝛼) 𝑆𝐸(𝛽)
𝑆𝐸 ( ) ∼ √( ) +( )
𝛽 𝛽 𝛼 𝛽

6. Retour aux paramètres de Michaelis-Menten

# Valeurs moyennes des paramètres
vmax <- 1/estim["(Intercept)"]
KM <- estim["x"]*vmax
[Link](vmax, KM, [Link]="Parameters")

vmax KM
Parameters 0.585323 0.4406275
# Erreurs standards
[Link] <- SE["(Intercept)"]/estim["(Intercept)"]^2
tmp <- sqrt((SE["(Intercept)"]/estim["(Intercept)"])^2+(SE["x"]/estim["x"])^2)
[Link] <- (estim["x"]/estim["(Intercept)"])*tmp
[Link]([Link], [Link], [Link]="[Link]")

[Link] [Link]
[Link] 0.1038999 0.09061058
# Calcul de la précision = erreur relative
[Link] <- 100*qt(0.975, n-2)*[Link]/vmax
[Link] <- 100*qt(0.975, n-2)*[Link]/KM
[Link]([Link], [Link], [Link]="ER(%)")
68 CHAPTER 4. VITESSE D’UNE RÉACTION ENZYMATIQUE

[Link] [Link]
ER(%) 49.2843 57.09478

• Qu’en pensez-vous ?

Vous voilà prêts maintenant à apprendre les instructions R pour faire de la

régression non linéaire.

4.4.2 Modèle non linéaire gaussien avec R

Nous allons donc maintenant revenir au modèle de Michaelis-Menten dans sa
formulation initiale :
𝑣𝑚𝑎𝑥 [𝑆]
𝑣([𝑆]) =
𝐾𝑀 + [𝑆]

Nous garderons le jeu de données n°1 pour illustrer les différentes instructions.
Nous utiliserons la librairie nlstools du logiciel R qui met à disposition tout un
panel d’outils de visualisation et d’analyse des résultats (Baty et al., 2015).

1. Prévisualisation du modèle

Comme évoqué plus haut (§ 4.2.3), les algorithmes de minimisation du critère

des moindres carrés, nécessitent de donner aux paramètres des valeurs initiales.
On peut faire cela “à l’oeil”, ou bien se référer à la signification biologique et/ou
géométrique des paramètres (Fig. 5). On peut donc s’aider des données elles-
mêmes pour trouver des valeurs initiales des paramètres (Fig. 11). Ainsi, pour
le jeu de données 1, on a 𝑣𝑚𝑎𝑥 ≈ 0.5 et 𝐾𝑀 ≈ 0.3.

Vérifions avec R qu’on n’est pas loin des valeurs optimales.

library(nlstools)

# ATTENTION : choisir le chemin d'accès vers votre propre fichier

dataMM1 <- [Link]("[Link]", header=TRUE)

n <- nrow(dataMM1) # nombre d'observations

# Implémentation du modèle
# ATTENTION : le nom des variables doit correspondre exactement
# aux en-têtes de colonne du jeu de données. Ici S et V.
modeleMM <- [Link]("V ~ vmax * S/(KM + S)")
# valeurs initiales des paramètres
valinit <- list(vmax = 0.5, KM = 0.3)
par(mfrow = c(1, 1))
par(mar = c(4, 4, 0.2, 0.2))
preview(formula = modeleMM, data = dataMM1, start = valinit)
4.4. RETOUR À LA RÉGRESSION GAUSSIENNE 69

0.50
0.45

+
0.40

+
0.35
Predicted

+
0.30
0.25

+
0.20
0.15

0.2 0.4 0.6 0.8 1.0 1.2 1.4

La valeur RSS: 0.00895 fournie par la fonction preview() après le graphe corre-
spond à la valeur de la 𝑆𝐶𝐸 pour les valeurs initiales choisies (valinit).

2. Ajustement par régression non linéaire

ajusMM <- nls(formula = modeleMM, data = dataMM1, start = valinit)
# Pour visualiser le graphe de l'ajustement :
plotfit(ajusMM, smooth = TRUE)
70 CHAPTER 4. VITESSE D’UNE RÉACTION ENZYMATIQUE

3. Premiers résultats On récupère tout un ensemble de résultats liés à

l’ajustement :
• Les paramètres estimés : valeurs moyennes (Estimate) et écart-types (Std.
Error) ; on obtient la même information avec l’instruction coef(ajusMM).
• L’estimation de 𝜎 (Residual standard error) ;
• La valeur minimum de la 𝑆𝐶𝐸 (Residual sum of squares) ;
Vérifiez numériquement la relation entre 𝜎 et 𝑆𝐶𝐸.
• Et d’autres résultats dont nous ne parlerons pas.
overview(ajusMM)

------
Formula: V ~ vmax * S/(KM + S)

Parameters:
Estimate Std. Error t value Pr(>|t|)
vmax 0.69040 0.03682 18.75 4.77e-05 ***
KM 0.59655 0.06826 8.74 0.000944 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.01356 on 4 degrees of freedom

4.4. RETOUR À LA RÉGRESSION GAUSSIENNE 71

Number of iterations to convergence: 5

Achieved convergence tolerance: 1.764e-06

------
Residual sum of squares: 0.000736

------
t-based confidence interval:
2.5% 97.5%
vmax 0.5881579 0.7926407
KM 0.4070356 0.7860596

------
Correlation matrix:
vmax KM
vmax 1.0000000 0.9450641
KM 0.9450641 1.0000000
4. Analyse des résidus - graphiques
residus <- nlsResiduals(ajusMM)
par(mar = c(4, 4, 3, 0.2))
plot(residus)
Residuals Standardized Residuals
Standardized residuals

2
Residuals

0.005

0
−2
−0.015

0.2 0.3 0.4 0.5 0.2 0.3 0.4 0.5

Fitted values Fitted values

Normal Q−Q Plot of
Autocorrelation
Standardized Residuals
Sample Quantiles

1.0
Residuals i+1

0.000

0.0
−0.015

−1.0

−0.015 −0.005 0.005 0.015 −1.0 −0.5 0.0 0.5 1.0

Residuals i Theoretical Quantiles

• Qu’en pensez-vous ?
Les résidus standardisés sont les résidus divisés par l’écart-type 𝜎. On constate
72 CHAPTER 4. VITESSE D’UNE RÉACTION ENZYMATIQUE

ici une certaine auto-corrélation

4.4.3 Comparaison des deux approches

Nous allons maintenant comparer les valeurs des paramètres estimés par linéari-
sation de Lineweaver & Burk et par régression non linéaire. Pour faire cette
comparaison, nous pouvons par exemple nous appuyer sur l’erreur relative des
estimations et sur des graphiques montrant simultanément les deux ajustements.
Nous aurons besoin des objets créés précédemment sous R, en particulier m et
ajusMM.
Pour comparer les erreurs relatives, nous devons travailler directement sur les
paramètres 𝑣𝑚𝑎𝑥 et 𝐾𝑀 .
# Pour la linéarisation
# On a déjà calculé [Link] et [Link]

# Pour la régression non linéaire

[Link] <- summary(ajusMM)$coefficients[,"Estimate"]
[Link] <- summary(ajusMM)$coefficients[,"Std. Error"]
[Link] <- 100*qt(0.975, n-2)*[Link]/[Link]

# Présentation de la comparaison sous forme d'une matrice

comp <- matrix(c([Link], [Link], [Link]), ncol=2, nrow=2)
colnames(comp) <- c("L&B", "MM")
rownames(comp) <- c("vmax", "KM")
print(comp)

L&B MM
vmax 49.28430 14.80903
KM 57.09478 31.76813
• Que pouvez-vous en conclure ?
Pour les graphiques, on peut choisir soit le plan ([𝑆], 𝑣) du modèle de Michaelis-
Menten, soit le plan (1/[𝑆], 1/𝑣) de la linéarisation de Lineweaver & Burk (voir
figure ci-dessous).
par(mfrow = c(1,2), mar = c(4,4,2,1))

# Dans le plan ([S], v)

plot(dataMM1$S, dataMM1$V, xlab="[S]", ylab="v", las=1, pch=19, main="Représentation di
curve([Link]["vmax"]*x/([Link]["KM"]+x), from=0, to=max(dataMM1$S), add=TRU
curve(vmax*x/(KM+x), from=0, to=max(dataMM1$S), add=TRUE, lty=2)
legend("bottomright", legend=c("MM", "L&B"), lty=c(1,2), bty="n", cex=0.75)

# Dans le plan (1/[S], 1/v)

plot(1/dataMM1$S, 1/dataMM1$V, xlab="1/[S]", ylab="1/v", las=1, pch=19, main="Représent
4.5. A VOUS DE JOUER ! 73

abline(a=1/coef(ajusMM)["vmax"], b=coef(ajusMM)["KM"]/coef(ajusMM)["vmax"])
abline(m, lty=2)
legend("bottomright", legend=c("MM", "L&B"), lty=c(1,2), bty="n", cex=0.75)
Représentation directe Représentation L&B

0.5 7

6
0.4
5

0.3 4

1/v
v

3
0.2

2
0.1
1
MM MM
0.0 L&B 0 L&B

0.0 0.5 1.0 1.5 0 2 4 6

[S] 1/[S]

• Que pouvez-vous en conclure ?

On voit sur ces deux graphes que la transformation de L&B donne un plus mau-
vais ajustement que la régression directe avec le modèle de MM, en particulier
pour les fortes valeurs de concentrations en substrat.

4.5 A vous de jouer !

Télécharger les jeux de données 2 et 3.
Pour chaque jeu de données, vous ajusterez les trois transformations linéaires
ainsi que la version non linéaire du modèle de Michaelis-Menten.
Pour chaque ajustement :
• vous représenterez graphiquement le modèle superposé aux données ;
• vous donnerez l’estimation moyenne des paramètres avec l’intervalle de
confiance à 95 % ;
• vous jugerez de la corrélation entre les paramètres ;
• vous examinerez les résidus ;
74 CHAPTER 4. VITESSE D’UNE RÉACTION ENZYMATIQUE

• vous vérifierez l’hypothèse de normalité des résidus.

Sur la base de vos résultats, vous discuterez du bien fondé des transformations
linéaires.
Vous veillerez également à toujours comparer les résultats de l’ajustement avec
le modèle de Michaelis-Menten, aux résultats des ajustements avec les différentes
transformations linéaires.
Eléments de corrections pour le modèle de Michaelis-Menten.
• Atkins G.L. et Nimmo I.A. 1975. A comparison of seven methods for
fitting the Michaelis-Menten equation. Biochemical Journal, 149:775-777.
• Baranyi J. et Terry R.A. 1995. Mathematics of predictive food microbiol-
ogy. International Journal of Food Microbiology, 26(2):199-218.
• Baty F., Ritz C., Charles S., Brutsche M., Flandrois J.-P., Delignette-
Muller M.L. 2015. The R package nlstools : a toolbox for nonlinear re-
gression. Journal of Statistical Software, 66:1–21.
• Bernard O. 2004. La modélisation des systèmes biologiques : Aller-retours
le long des fleuves qui circulent entre l’océan du réel et le lac des modèles.
Habilitation à diriger des recherches, Université de Nice-Sophia-Antipolis.
105p.
• Buchanan R.E. 1918. Life phases in a bacterial culture. The Journal of
Infectious Diseases, 23:109–123.
• Eadie G.S. 1942. The inhibition of cholinesterase by physostigmine and
prostigmine. Journal of Biological Chemistry, 146:85-93.
• Edelstein-Keshet L. 2005. Mathematical Models in Biology. Classics in
Applied Mathematics, SIAM, Philadelphia, 586p.
• Galton F. 1890. Kinship and Correlation. The North American Review,
150:419–431.
• Hofstee B.H.J. 1952. On the evaluation of the constants 𝑉𝑚 and 𝐾𝑀 in
enzyme reactions. Science, 116:329-331.
• Lineweaver H. et Burk D. 1934. The determination of enzyme dissociation
constants. Journal of the American Chemical Society, 56:658-666.
• Lobry J.R. 1991. Ré-évaluation du modèle de croissance de Monod. Effet
des antibiotiques sur l’énergie de maintenance. Thèse de doctorat “Ecolo-
gie, Environnement”. Université Claude Bernard - Lyon I, 178p..
• Malthus T. 1798. The Principle of Population. London. ©1998, Electronic
Scholarly Publishing Project.
• Michaelis L. et Menten M.L. 1913. Die Kinetik der Invertinwirkung. Bio-
chemische Zeitschrift, 49:333–369.
4.5. A VOUS DE JOUER ! 75

• Monod J. 1935. Le taux de croissance en fonction de la concentration de

l’alimnent dans une population de Glaucoma piriformis en culture pure.
Comptes-Rendus l’Académie des Sciences, 201:1513–1515.
• Monod J. 1942. Recherches sur la croissance des cultures bactériennes,
Herman éd., Paris.
• Monod J. 1949. The Growth of Bacterial Cultures. Annual Reviews of
Microbiology, 3:371–394.
• Monod J, Tessier G. 1936. La concentration de l’aliment, facteur quan-
titatif de l’accroissement des populations d’Infusoires. Comptes-Rendus
l’Académie des Sciences, 202:162–164.
• Novick A., Szilard L. 1950. Description of the Chemostat. Science,
112:715–716.
• Pavé A. 2012. Modélisation des systèmes vivants. Lavoisier, Paris. 633p.
• Pletcher SD. 1999. Model fitting and hypothesis testing for age-specific
mortality data. Journal of Evolutionnary Biology, 12:430–439.
• Ritchie RJ, Prvan T. 1996. Current statistical methods for estimating
the Km and Vmax of Michaelis-Menten kinetics. Biochemical Education,
24:196–206.
• Verhulst P. 1845. Recherches mathématiques sur la loi d’accroissement de
la population. Nouveaux mémoires de l’Académie Royale des Sciences et
Belles-Lettres de Bruxelles. T. XVIII, 25p.
• Wilkinson G.N. 1961. Statistical estimations in enzyme kinetics. Bio-
chemical Journal, 80:324-332.
• Zwietering M.H., Jongenburger I., Rombouts F.M., Van’t Riet K. 1990.
Modeling of the bacterial growth curve. Applied Environmental Microbi-
ology, 56:1875–1881.