phd-thesis-fr/partie-otb.tex

\chapter{Travelling Bacteria}\label{chap:partie-otb}
\minitoc
\begin{figure}[h]
  \centerline{%
    \includegraphics[width=19cm]{figures/otb2-display}}
  \caption{Fenêtre de rendu du logiciel \tb{} lors d'une simulation comportant
    un grand nombre de bactéries. Les cadres jaunes indiquent l'espace
    actuellement occupé par les bactéries et leur support}
  \label{fig:otb2-display}
\end{figure}

\clearpage
\section*{Préliminaires}
Voici quelques notations utilisées dans ce chapitre:
\begin{itemize}
\item $(G,+)$ est le groupe engendré par l'ensemble $G$ et la loi $+ : G \times
  G \rightarrow G$.
\item $\langle G \rangle$ est le groupe libre engendré par l'ensemble des
  générateurs $G$.
\item $\mathbb{N}_+$ est l'ensemble des entiers naturels privé de $0$.
\item $\lfloor x \rfloor = \text{max} \{ m \in \mathbb{Z} \,|\, m \leq x \}$ est
  la fonction partie entière de $x \in \mathbb{R}$.
\item $\orr{u}$ est un vecteur,
\item $\|\orr{u}\|$ est la norme euclidienne du vecteur $\orr{u}$,
\item $k \cdot\orr{u}$ est le produit scalaire entre un scalaire $k$ et un
  vecteur \orr{u},
\item $\orr{u}\cdot\orr{v}$ est le produit scalaire entre deux vecteurs \orr{u}
  et \orr{v},
\item $\orr{u}\wedge\orr{v}$ est le produit vectoriel entre deux vecteurs
  \orr{u} et \orr{v}.
\item $\left[ x \right]^{x_{\mathit{max}}}_{x_{\mathit{min}}}
  = \left\{ \begin{array}{ll}
              x_{\mathit{max}} & \text{si}\; x > \mathit{max}\\
              x_{\mathit{min}} & \text{si}\; x < \mathit{min}\\
              x & \text{sinon}
            \end{array}
          \right.$ est la fonction \emph{clamp}.
\end{itemize}

\section{Simuler la morphogenèse des populations}
En tant que science expérimentale, la recherche en biologie cellulaire et
moléculaire nécessite d'important moyens matériels afin de tester \emph{in
vivo} les hypothèses formulées par les chercheurs. Notamment, la préparation
et l'exécution d'une expérience est…
\begin{itemize}
\item consommatrice en temps. Il faut parfois attendre que les cultures d'un
  organisme particulier atteignent l'état désiré\footnote{cf. le projet
    \synbiotic};
\item consommatrice en ressources. Maintenir les conditions de pression,
  d'humidité, de température, d'acidité, d'oxygénation, etc. nécessite du
  matériel onéreux;
\end{itemize}

D'un côté, de par les moyens engagés, il est difficile d'envisager de tester des
hypothèses très prospectives ou bien de se faire une idée du résultat probable
d'une expérience quand beaucoup d'entités interagissent. D'un autre, un plan
d'expérience comprend idéalement un nombre minimum nécessaire d'expérimentations
pour valider ou infirmer une hypothèse du modèle, ce nombre pouvant néanmoins
rester élevé au regard des moyens à disposition. La simulation informatique
permet, à partir d'un modèle donné, de confirmer ou de rejeter rapidement une
hypothèse, ce qui peut aider à choisir quelles expérimentations effectuer. Dans
ce cas c'est un outil d'aide à la décision. Si la simulation est suffisamment
précise, comme dans la cas de la simulation de la motilité de quelques bactéries
\Ecoli{} \cite{bray_chemotactic_2007}, elle pourrait même se substituer en
partie à l'expérimentation \emph{in vivo}; on parlera dans ce cas
d'expérimentation \emph{in silico}.

Nous nous intéressons au problème de la \emph{morphogénèse} c'est à dire à
l'ensemble des lois qui déterminent les étapes de la forme, de la structure des
tissus, des organes et des organismes. Les organismes multicellulaires sont le
fruit d'un grand nombre de divisions cellulaires successives et qui donnent lieu
à l'apparition d'une forme et de fonctionnalités caractéristiques et
reproductibles. Comprendre la morphogénèse est un enjeu fondamental pour la
biologie synthétique (voir à ce propos la section~\ref{sec:synbiotic}).

Pour répondre à ce problème et dans le cadre du projet \synbiotic, nous avons
développé \tb{}, un simulateur des interactions d'une population de bactéries
\ecoli{} au cours de sa croissance. \tb{} permet de suivre une cellule ou une
sous partie de cette population au cours du temps de manière similaire aux
conditions expérimentales et ainsi de tester des hypothèses à l'échelle de la
population dans le domaine de la biologie synthétique.

\subsection{Propriétés attendues}
Dans le but de fournir aux biologistes le moyen d'expérimenter \emph{in silico}
la validité d'un modèle concernant l'évolution d'une population de bactéries
\ecoli{}, nous avons dressé une liste de propriétés nécessaires.
\begin{description}
\item[Sorties] Le simulateur permet de recueillir des résultats qualitatifs via
  un affichage graphique, reprenant la métaphore de la boîte de Petri qui
  permettra d'identifier la forme de la population en croissance. Il est
  possible de marquer une sous-partie de la population pour la suivre plus
  facilement. Le simulateur permet de recueillir des résultats quantitatifs, tel
  que le nombre et la position des bactéries dans un fichier texte pour
  traitement ultérieur, par exemple pour la reconnaissance automatique de forme,
  ou des calculs statistiques;
\item[Entrées] Les conditions initiales de la simulation sont décrites par un
  fichier de configuration dans un langage dédié; au delà des conditions
  initiales de l'expérience, il permettra au moins de fixer quelques paramètres
  du rendu graphique ou textuel, comme la couleur de représentation des
  morphogènes et des bactéries, et de sélectionner quelles informations
  quantitatives seront écrites pendant la simulation;
\item[Précision] Le simulateur reprend aussi fidèlement que possible les
  connaissances actuelles sur l'aspect mécanique (croissance et collision des
  bactéries) et sur l'aspect chimique (réaction et diffusion de morphogènes).
  Une étape de calibrage sera donc nécessaire par rapport à des comptes rendus
  d'expériences \emph{in vivo};
\item[Efficacité] Le simulateur est suffisamment efficace pour être utilisable
  sans matériel spécifique; par exemple, passer d'une bactérie à une population
  de l'ordre de $10^4$ bactéries, soit un temps réel d'environ de
  \SI{4}{\hour}\SI{30}{\minute}, nécessitera un temps de simulation
  correspondant de l'ordre de l'heure sur un ordinateur portable contemporain,
  moins encore sur une machine spécialisée;
\item[Autonomie] Les résultats sont accessibles sans la nécessité de faire
  fonctionner le simulateur sur une grille de calcul, un tel matériel ne servira
  qu'à réduire le temps de simulation afin d'effectuer un plan d'expérience plus
  rapidement;
\item[Scalabilité] Malgré les points précédents, \tb{} passe à l'échelle et peut
  s'adapter à une machine de calcul dédiée ou à une grille de calcul distribuée;
  plusieurs instances du simulateur peuvent être lancées de manière concurrente;
\item[Adaptabilité] Il est possible d'ajouter à \tb{} de nouvelles
  fonctionnalités avec peu d'effort: l'architecture du simulateur est modulaire
  et vient avec une documentation claire;
\item[Stabilité] Comme une simulation complexe peut prendre plusieurs heures, il
  est important que le simulateur ne produise pas d'erreur en cours de
  fonctionnement, autrement dit \tb{} est stable; d'une part le langage utilisé
  apportera des garanties de part des vérifications du système de type, d'autre
  part des tests fonctionnels permettront de valider le fonctionnement correct
  de chaque partie du simulateur;
\item[Ouverture] Le simulateur est un programme informatique ouvert, son code
  source peut-être téléchargé en ligne, compilé, amélioré et redistribué,
  modifié ou non;
\end{description}

\subsection{\ocaml} % Adaptabilité, stabilité
Pour répondre aux critères \emph{d'adaptabilité} et de \emph{stabilité},
nous avons choisi d'implémenter \otb avec le langage \ocaml, un langage de
programmation fonctionnel, statiquement typé, modulaire et interfaçable avec
le langage C. \ocaml est à la fois un projet académique avec de solides
bases théoriques, notamment utilisé comme langage d'implémentation de
l'assistant de preuve Coq\footnote{cf. la page «à propos» du projet Coq
\url{https://coq.inria.fr/about-coq}}, et aussi un langage industriel ayant à
son actif de belles réussites comme l'analyseur statique \textsc{Astrée} en
usage chez Airbus.

\paragraph{Langage de programmation fonctionnel}
La programmation fonctionnelle est issue du $\lambda$-calcul, un modèle de
calcul mathématique reposant sur des variables, des applications et des
abstractions regroupées sous le nom de \emph{$\lambda$-termes}.
Dans le modèle du $\lambda$-calcul, la règle centrale pour le calcul des
$\lambda$-termes est la \emph{substitution} dont les deux principaux
représentants sont l'$\alpha$-conversion et la $\beta$-réduction. Le calcul de
nouveaux $\lambda$-termes s'effectue par un système de réécriture. C'est ce
modèle de calcul qui a donné lieu au paradigme de la programmation
fonctionnelle. Par extension, les langages de programmation fonctionnelle purs
ont la particularité de considérer la programmation d'un algorithme comme une
composition de fonctions s'appelant les unes les autres plutôt que comme une
suite d'instructions qui agissent par effet de bord, c'est-à-dire en modifiant
l'état de la mémoire. Ils se dissocient de la vision de la machine à état comme
support de programmation en adoptant un style plus mathématique dans leur
écriture, proche de la réécriture des $\lambda$-termes. \ocaml dispose ainsi
\begin{itemize}
\item de structures de contrôle comme le filtrage par motif,
\item de clôtures, qui sont des définitions de fonction accompagnées de leur
  environnement lexical, analogues des $\lambda$-termes et
\item permet de définir des fonctions d'ordre supérieur, c'est-à-dire aux
  fonctions pouvant avoir des fonctions en argument et pouvant retourner des
  fonctions comme résultat.
\end{itemize}
Un paradigme de programmation influe sur la syntaxe et la sémantique du langage,
et comme \ocaml est un langage multiparadigme, empruntant aux paradigmes
de la programmation fonctionnelle, de la programmation impérative et de la
programmation orientée objet, il hérite de leurs aspects. \ocaml dispose d'un
ramasse-miette efficace qui gère la désallocation de la mémoire quand celle-ci
n'est plus utilisée sans l'intervention de l'utilisateur. Ceci fait de \ocaml
un langage de programmation de haut niveau. Comme nous le présentons ci-après,
\ocaml est de plus un langage fonctionnel fortement et statiquement typé.

\paragraph{Typage fort et typage statique}
Les types ont d'abord été introduit en logique combinatoire et dans le
$\lambda$-calcul dans les années 1930. Ils ont depuis eu un fort impact en
informatique: ils sont utilisés comme une abstraction d'un programme et
permettent d'en vérifier certaines propriétés, comme par exemple garantir que le
résultat d'une fonction appartient à un ensemble d'éléments connus. Dans le
langage C, un langage bas niveau proche de l'expression d'un processeur, les
types décrivent la taille de la représentation machine et les règles de
conversion d'une valeur du langage. En \ocaml, les types sont plutôt utilisés
pour modéliser les données du programme. \ocaml permet de déclarer des types
algébriques récursifs permettant par exemple de décrire la structure de listes
\begin{OCAMLcode}
type 'a list =
    Nil
  | Cons of 'a * 'a list
\end{OCAMLcode}
ou d'arbres binaires
\begin{OCAMLcode}
type 'a tree =
    Leaf of 'a
  | Node of 'a tree * 'a tree
\end{OCAMLcode}

\ocaml présente un typage fort. On distingue un typage fort d'un typage faible
par la manière dont s'effectue la transformation d'un type en un autre: plus
elle est explicite plus c'est un typage fort, moins elle explicite moins c'est
un typage faible. Autrement dit, plus un langage est fortement typé, plus il est
difficile de transformer un type en un autre. Dans le langage C, un opérateur
appelé \emph{cast} permet aisément de transformer presque n'importe quel type en
un autre. En \ocaml, cette transformation n'est pas possible dans le cadre de
fonctionnement normal du langage, donc \ocaml est plus fortement typé que C et
le système de type d'\ocaml est par conséquent plus riche. Le code source d'un
programme \ocaml est typé statiquement, c'est à dire que toutes les
annotations de type sont connues avant la compilation. À l'opposé, dans le cas
d'un typage dynamique, certaines informations de type peuvent ne pas être
connues avant la compilation, peuvent rester après la compilation et peuvent
évoluer durant l'exécution du programme.

\ocaml utilise un typage statique et fort et les annotations de type sont
utilisées à la compilation pour permettre au compilateur de détecter les erreurs
de types avant l'exécution du programme, ce qui apporte quelques garanties
concernant le comportement. Par exemple, les types nous permettent d'apporter
des informations supplémentaires que le compilateur peut prendre en compte.
Déclarer les types
\begin{OCAMLcode}
  type nom = Nom of string
  type prenom = Prenom of string
\end{OCAMLcode}
permet au compilateur de distinguer deux chaînes de caractères identiques, mais
dont le sens est différent, ce qui est représenté par deux types différents.

\paragraph{Modularité native}
Lors de la rédaction d'un programme, il est vite nécessaire de «compartimenter»
le code source, d'abord parce que certaines fonctionnalité sont génériques et
peuvent être utilisées plusieurs fois dans des contextes différents, mais aussi
pour des raisons de relecture et de segmentation de l'espace de noms. La
\emph{modularité} consiste à rassembler les fonctionnalités d'un programme dans
des modules \emph{indépendants} et \emph{interchangeables}. En \ocaml, la
modularité est une caractéristique native du langage dans la gestion du code
source et repose sur trois constructions syntaxiques: les \emph{structures}, les
\emph{interfaces} et les \emph{modules}. Les structures contiennent
l'implémentation des modules, les interfaces décrivent les valeurs qui en sont
accessibles — les valeurs dont l'implémentation n'est pas exposée sont des
valeurs abstraites, et celles qui n'apparaissent pas du tout dans
l'implémentation du module sont inaccessibles. Un module peut avoir plusieurs
interfaces — du moment qu'elles sont toutes compatibles avec les types de
l'implémentation — et plusieurs modules peuvent vérifier une même interface.
Enfin, \ocaml dispose de modules paramétriques nommés \emph{foncteurs} qui
sont des structures paramétrées par d'autres structures.

% http://caml.inria.fr/pub/docs/manual-ocaml/intfc.html
\paragraph{Interface avec le langage C}
Un aspect important pour les besoins du cahier des charges (voir
\ref{ss:req-parallel}) est de pouvoir communiquer avec des objets compilés C via
une interface dédiée, typiquement des bibliothèques écrites et compilés dans ce
langage. \ocaml propose une interface claire avec des bibliothèques
logicielles écrites en C, permettant d'empaqueter des valeurs dans des types
\ocaml conservant les avantages d'un typage fort, de dépaqueter des valeurs
typées de \ocaml vers un type C arbitraire et surtout de profiter de la
gestion automatique de la mémoire par le ramasse-miette de \ocaml.

\subsection{Calcul parallèle} % Efficacité et Autonomie
\label{ss:req-parallel}
Notre simulateur calcule l'évolution d'un grand nombre d'objet de même nature
simultanément. Pour répondre aux critères \emph{d'efficacité} et
\emph{d'autonomie} nous avons le choix d'utiliser un GPU, le processeur dédié
des cartes graphiques, ou bien d'utiliser un cluster: une grappe de serveurs
interconnectés sur le réseau local. En prenant le nombre d'opérations en virgule
flottante par seconde, notée \si{FLOPS}, comme référence de la vitesse d'un
ordinateur, utiliser le parallélisme embarqué dans les GPU semble être la bonne
manière d'atteindre des vitesses de calcul élevées. Nous nous appuyons sur le
fait que plus de la moitié des 500 supercalculateurs les plus rapides au monde
sont construit sur NVidia
Kepler\footnote{\url{https://www.top500.org/statistics/list/}}, un GPU. Les GPU
grand public sont accessibles en terme de prix, leur consommation est
suffisamment faible pour fonctionner sur batterie dans un PC portable et les
plus puissants atteignent les \SI{10000}{\giga FLOPS}\footnote{%
  \url{https://www.nvidia.com/en-us/geforce/products/10series/titan-x-pascal}\\
  Par exemple, la carte graphique NVidia Titan X possède \num{3584} cœurs
  de calculs, cadencés à \SI{1417}{\mega\hertz} pour un prix avoisinant les
\num{1200} euros}.Enfin,
le parallélisme des GPU semble être plus efficace énergétiquement par rapport
aux CPU en comparant leur \si{\giga FLOPS} par Watt~\cite{dong_step_2014}.

%parallélisme SIMD OpenCL
\paragraph{SIMD}
La taxinomie de Flynn~\cite{flynn_computer_1972} classifie différentes
architectures d'ordinateur et comporte quatre catégories selon le type
d'organisation du flux de données et du flux d'instructions. Le modèle de
programmation le plus répandu repose sur l'architecture de Von Neumann, ce qui
correspond à la catégorie \emph{Single Instruction Single Data} (SISD). Les
ordinateurs de cette catégorie n'exploitent aucun parallélisme et sont purement
séquentiels, tant au niveau des instructions qu'au niveau de la mémoire. Ce
modèle de programmation ne correspond pas à la réalité des processeurs
contemporains qui exploitent le parallélisme des données (plusieurs banques de
mémoire RAM) ou d'instructions grâce à plusieurs processeurs à plusieurs cœurs.
À l'extrême opposé de la classification se trouvent les architectures
\emph{Multiple Instruction on Multiple Data} (MIMD), puis les architectures
\emph{Multiple Instructions on Single Data} (MISD). Dans le cas des MISD, des
problèmes de concurrence à la lecture−écriture apparaissent et il est d'usage
d'avoir recours à des barrières de synchronisation, des sémaphores ou des Mutex.
SPMD, \emph{Single Program on Multiple Data}, est la technique la plus courante
pour utiliser le parallélisme des ordinateurs multiprocesseurs contemporain, où
chaque processeur exécute un programme distinct et possède sa propre mémoire, en
plus de partager une mémoire globale, la RAM. Finalement, dans les ordinateurs
contemporains, il existe, en complément d'un CPU et d'une mémoire partagée, des
cartes de calcul dédiées, comme des cartes graphiques, spécialisées dans le
calcul vectoriel, capable d'effectuer des opérations rapides sur des grandes
matrices. L'architecture d'une carte graphique, munie d'un GPU et d'une mémoire
RAM entre dans la catégorie \emph{Single Instruction on Multiple Data} (SIMD).
Ces processeurs sont spécialisés dans le traitement en parallèle d'un grand
nombre de données. Depuis la taxinomie de Flynn, d'autres types de Nous avons
choisi de tirer parti du parallélisme disponible sur l'ordinateur où est exécuté
\tb{}, qu'il provienne de plusieurs CPU ou d'une carte graphique.

\paragraph{\opencl}
Afin d'accéder aux périphériques présent sur la plateforme de simulation et de
les programmer, nous utilisons \opencl (Open Computing Language), un framework
pour l'écriture de programmes s'exécutant sur des plateformes hétérogènes:
microprocesseurs (CPU), processeurs graphique (GPU), processeurs de signal
numérique (DSP), réseau de portes programmables in situ (FPGA) ou autre carte
accélératrice. \opencl spécifie un langage de programmation, inspiré de C99,
pour la programmation de ces périphériques ainsi que des interfaces de
programmation (API) pour le contrôle de la plateforme et l'exécution des
programmes sur les périphériques. \opencl fournit une interface standard pour
la programmation parallèle utilisant le parallélisme de tâche et le parallélisme
de données. Il existe d'autres framework pour la programmation des GPU, comme
CUDA\footnote{\url{http://www.nvidia.com/object/cuda_home_new.html}}, cependant
l'avantage de \opencl est qu'il est le seul langage pour le GPGPU (General
Purpose computing on Graphics Processing Unit) à ne pas être restreint à une
seule plateforme. En effet, il fonctionne sur une large gamme de matériel et, au
contraire de CUDA, ne se limite pas à un seul constructeur, ni à un seul type de
périphérique. \opencl est un standard ouvert maintenu par Khronos
Group\footnote{\url{https://www.khronos.org/about}}, un consortium technologique
à but non lucratif. De nombreux constructeurs fournissent une implémentation
officielle, comme AMD, Intel et Nvidia les principaux constructeurs de CPU et de
cartes graphiques grand public.

\paragraph{\opengl} % Sorties
Pour répondre au point «Sortie», nous avons choisi d'utiliser \opengl (Open
Graphics Library), une API multiplateforme, multilangage utilisée pour
programmer un GPU spécifiquement pour un rendu graphique. Dans le modèle de
programmation de \opengl, un ensemble de vecteurs dans l'espace en trois
dimensions sont successivement transformés pour obtenir une projection en deux
dimension rastérisée. Le programmeur écrit des programme de traitement, des
shader, qui viennent se brancher à différents points de la chaîne de traitement
graphique. À l'inverse de \opencl, \opengl est optimisé pour le calcul vectoriel
et n'est pas approprié pour effectuer des calculs généraux sur les GPU. Cette
bibliothèque logicielle nous permet de proposer une ou plusieurs fenêtre
de rendu en temps réel pour présenter des informations qualitatives sur la
progression d'une simulation.

\subsection{Modélisation d'une cellule}
\begin{figure}[t]
  \centering
  \includegraphics[height=6cm]{figures/EscherichiaColi}\hfill
  \includegraphics[page=1,height=6cm]{figures/bacteria}
  \caption{Une image d'une microscopie électronique à balayage d'une culture de
    \ecoli\ (à gauche) et une représentation simplifiée en deux dimensions d'une
    bactérie \ecoli\ (à droite)}
  \label{fig:ecoli}
\end{figure}

\Ecoli\ est une bactérie bacillaire, c'est à dire en forme de bâtonnet, du genre
\emph{Escherichia}. C'est une bactérie vivant habituellement dans l'intestin des
animaux à sang chaud, dont elle constitue 0,1\% de la flore intestinale.
Découverte en 1885 par Theodor Escherich et étudiée intensivement depuis 60 ans,
c'est un organisme très utilisé dans le domaine de l'ingénierie biologique et de
la microbiologie industrielle, notamment pour sa disponibilité et sa facilité de
culture. La souche K-12, adaptée à l'étude en laboratoire, est également un des
premiers organisme dont le génome a été séquencé \cite{blattner_complete_1997}.

Elle est à l'origine du domaine de la biotechnologie en étant choisie comme
premier organisme génétiquement modifié contenant de l'ADN recombiné, c'est à
dire de l'ADN dont la séquence des nucléotides à été altérée pour modifier
l'expression des gènes de la bactérie. Elle est ensuite utilisée comme hôte pour
la production de protéines hétérologues (provenant d'un autre organisme),
notamment pour la production industrielle d'insuline humaine
\cite{goeddel_expression_1979} mais aussi de vaccins. \ecoli{} est utilisée en
bioremédiation pour la neutralisation de déchets toxique ou la dépollution des
sites contaminés et pour la production d'enzymes. Les efforts de recherche se
concentrent actuellement sur la synthèse par \ecoli{} de protéines plus
complexes comprenant par exemple des ponts bisulfure ou nécessitant une
modification post-traductionnelle pour être stables ou bien juste
fonctionnelles~\cite{han_escherichia_2006}.

\begin{table}[hb]
  \centering
  \begin{tabular}{llll}
    \hline
    \textbf{Paramètre}  & \textbf{Symbole}  & \textbf{Valeur}   & \textbf{Unité}\\
    \hline
    Longueur            & $L = 2 (l+r)$     & \num{2.5+-0.6}    & \si{\micro\meter} \\
    Largeur ou Diamètre & $d = 2r$          & \num{0.88+-0.09}  & \si{\micro\meter} \\
    Volume              & $V$               & \num{1.4}         & \si{\femto\liter} \\
    Vitesse rectiligne  & $v$               & \num{29+-6}       & \si{\micro\meter\per\second} \\
    Temps de génération & $G$               & 20 — 30           & \si{\minute}  \\
    \hline
  \end{tabular}
  \caption{Quelques ordres de grandeur associés aux bactéries \Ecoli{}
    \cite{bremer_modulation_1996,darnton_torque_2007}%
    \cite{britanica_online_encyclopedia_bacteria_2016}}
  \label{tab:ecoli-odg}
\end{table}

\paragraph{La chimiotaxie de \Ecoli}
Comme la plupart des bactéries bacillaires, \ecoli{} possède une ciliature
péritriche: elle est assortie de longs flagelles elliptiques, quatre en moyenne,
recouvrant sa surface qui sont utilisés comme moteur de son déplacement. Ces
flagelles sont fixés à leur base à un moteur rotatif réversible entraîné par un
flux de protons et tournent sur eux-mêmes à une fréquence de l'ordre de
\SI{100}{\hertz}. Suivant le sens de leur rotation, sens direct (CCW pour
Counter ClockWise) et sens indirect (CW pour ClockWise), ils sont à l'origine de
deux modes de déplacement:
\begin{itemize}
\item le mode «run», les moteurs tournent dans le sens direct et les flagelles
  se regroupent pour propulser la bactérie dans une direction rectiligne. Il
  dure en moyenne \SI{1}{\second};
\item le mode «tumble», au moins un des moteurs tourne dans le sens indirect.
  Les flagelles se désolidarisent ce qui a pour effet moyen de faire tourner la
  bactérie sur elle-même dans le sens horaire \cite{diluzio_escherichia_2005}.
  Ce mode dure en moyenne \SI{0,1}{\second}.
\end{itemize}
Ces deux modes de déplacement alternent successivement et, en fonction de
l'environnement de la bactérie, les amènent à migrer vers des environnements
favorables, riches en nutriments et à fuir les environnements toxiques et
délétères \cite{berg_chemotaxis_1972}.
Ce phénomène s'appelle la \emph{chimiotaxie}. Il résulte de la transmission
de signaux extracellulaires, détectés par les récepteurs de la membrane, aux
moteurs des flagelles qui contrôlent le comportement de la cellule.

D'un point de vue macroscopique, \ecoli{} effectue une marche aléatoire
biaisée vers une direction particulière. Elle détecte la variation
\emph{dans le temps} de la concentration d'un chimioeffecteur. Ainsi, lorsque
la bactérie se dirige vers un milieu plus favorable, les réorientations par
tumble deviennent plus rares et la bactérie peut remonter le gradient de la
concentration d'une molécule qui se diffuse dans l'environnement. Entre deux
run, \ecoli{} peut complètement changer de direction, avec un angle allant
jusqu'à 180°, le plus souvent entre 0° et 90° \cite{turner_real-time_2000}.
D'un point de vue microscopique, \ecoli{} peut détecter différents acides
aminés, sucres et dipeptides ainsi que l'acidité, la température et l'état
d'oxydoréduction.

Les récepteurs principaux, comme ceux de l'aspartate (Tar)
et de la sérine (Tsr), sont très nombreux et sont aux nombre de quelques
milliers par cellule. Les récepteurs secondaires, comme ceux spécifiques aux
dipeptides (Tap), ribose et galactose (Trg), et au potentiel d'oxydoréduction
(Aer), sont bien moins nombreux avec seulement une centaine de copies par
cellule. Ces récepteurs se regroupent en amas, ce qui induit des dynamiques
complexes dans la détection de chimioeffecteurs \cite{sourjik_receptor_2004},
impliquant soit une plus grande sensibilité ou au contraire une plus grande
diversité de détection.

Dans le réseau de régulation génétique, parmi les voies de signalisation
de la cellule, la voie de réponse chimiotaxique consiste donc en un ensemble
de protéines réceptrices transmembranaires sus-citées et les produits de six
gènes chimiotaxiques: CheR, CheB, CheW, CheA, CheY et CheZ. L'information sur
la liaison de ligants bénéfiques ou nocifs aux récepteurs est transportée
jusqu'au moteur des flagelles ce qui leur permet de changer le sens de rotation
du moteur. Il existe plusieurs modèles stochastiques capables de simuler
en détail la cascade d'évènements entre la liaison d'un chimioeffecteur
à un récepteur jusqu'au biais des moteurs \cite{shimizu_spatially_2003,
mello_quantitative_2003}.



\paragraph{Modèle de croissance}
\ecoli\ suit un modèle de croissance linéaire
\cite{zaritsky_growth_1982,kubitschek_cell_1990} par élongation et bien que son
diamètre décroisse avec l'élongation, cette variation reste négligeable quand
$L$ est suffisamment grande \cite{trueba_changes_1980,skarstad_cell_1983}. Nous
considérerons donc que seule sa longueur $L$ (voir \textsc{Fig.}~\ref{tab:ecoli-odg})
varie lors de la croissance d'une bactérie. Comme toutes les cellules
procaryotes, \ecoli{} se reproduit par scissiparité. Ce cycle de division
implique:
\begin{enumerate}
\item La réplication de la molécule d'ADN circulaire au sein de la cellule;
\item La migration après duplication des deux brins d'ADN à chacun des pôles de
  la cellule;
\item La croissance de la cellule pour laisser la place à un matériel génétique
  plus volumineux;
\item Le plan équatorial de la cellule se resserre et fini par scinder la
  membrane en deux, de sorte que chaque nouvelle cellule fille ait le même
  matériel génétique.
\end{enumerate}

L'équation \ref{eq:generation-time} présente le calcul du temps de génération
pour une population de bactéries dont la croissance s'effectue par division
binaire, où $t$ est le temps en minutes, $n$ est le nombre de générations, $B$
est le nombre de bactéries au début de la période de temps mesurée et $b$ est le
nombre de bactéries à la fin de la période de temps mesurée. Cela permet aux
biologistes de retrouver le temps moyen d'un cycle cellulaire pour la population
observée.

\begin{equation}
  \displaystyle
  G = \frac{t}{n} = \frac{t}{3,3 \times ln (b / B)}
  \label{eq:generation-time}
\end{equation}

\paragraph{Vieillissement et mort cellulaire}
Il semble que les bactéries se scindant en deux bactéries filles en tout point
identiques atteignent une sorte d'immortalité, au moins fonctionnelle. Néanmoins
quand elle n'est pas tuée par des antibiotiques, par manque de nutriments ou par
le système immunitaire d'organisme multicellulaires, même dans des conditions
idéales \ecoli{} vieillit~\cite{stewart_aging_2005}. Ce vieillissement résulte
d'une asymétrie entre les parois membranaires des pôles de la cellule. En effet,
un de ces pôle est celui de la cellule mère, alors que le second provient du
plan équatorial de scission de la cellule mère et a été fraîchement créé. En
partant d'une unique cellule mère, à chaque génération deux cellules possèdent
les parois membranaires des pôles de la cellule mère susmentionnée. Avec les
générations, le matériel cellulaire diffusant peu et ayant une longue demi-vie
s'accumule dans ces «vieux» pôles. Les cellules aux pôles plus âgés produisent
moins de cellules filles et ont une probabilité de décès spontané plus élevée.
On peut ainsi donner un âge à chaque pôle prenant en compte le nombre de
divisions qu'il a subit. Cela permet d'attribuer un âge à chaque cellule, qui
est définit comme l'âge du pôle le plus âgé. En conséquence, \otb inclus un
modèle de vieillissement des bactéries pouvant causer le décès spontané (nommé
apoptose) d'une cellule si sa lignée dépasse une certaine valeur, paramétrable
en début de simulation.

\subsection{Organisation logicielle du simulateur}
Pour décrire l'organisation logicielle du simulateur, il nous faut parler des
architectures de \opencl et de \opengl avant de présenter l'architecture de
\otb afin de comprendre comment ces différents composants communiquent pour
simuler l'évolution d'une population de bactéries.

\subsubsection{Architecture \opencl}\label{sec:archi-cl}
Nous présentons ici une version simplifiée de l'architecture de \opencl%
\footnote{\url{http://www.khronos.org/registry/cl/specs/opencl-1.1.pdf}}. Un
programme \opencl est construit en deux parties: le programme hôte et les
\emph{kernel}. Le programme hôte s'exécute normalement sur CPU. Dans notre cas,
c'est la partie du simulateur gérée par \ocaml. Les kernels sont des fonctions
écrites en C99 destinées à être exécutées sur les périphériques gérés par
\opencl, comme une carte graphique. Le modèle d'exécution d'\opencl définit
comment ces kernels sont instanciés sur les différentes unités de calcul
disponibles.

Lorsque l'hôte soumet un kernel à exécuter à \opencl, un espace d'indices est
définit et, pour chaque indice de cet espace, une instance du kernel est
exécutée. Cette instance est appelée un \emph{work-item}. Les work-items sont
regroupés en \emph{work-group}, une division par blocs de l'espace d'indices.
Les work-items ont accès à quatre types de mémoire: la mémoire globale, la
mémoire constante, la mémoire locale et la mémoire privée. Chacune de ces
mémoires est manipulable soit par l'hôte, soit par un work-item, soit par les
deux. Par exemple, la mémoire privée est réservée à un work-item, la mémoire
local est accessible à tous les work-items du même work-group, les mémoires
constantes et globales sont accessibles à l'hôte comme aux work-items.

L'hôte peut également demander des transferts de données entre sa mémoire
principale, la RAM, et les mémoires du device, dans les deux sens, afin de
fournir des données initiales aux work-items et de récupérer leur travail.

\subsubsection{Architecture \opengl}
Le framework \opengl définit une architecture client-serveur: le client
fonctionnant sur l'hôte, côté CPU, et le serveur s'exécutant sur la carte
graphique, côté GPU. La pièce maîtresse du serveur est le \emph{rendering
  pipeline} (chaîne de traitement du rendu), une succession de tâches permettant
de transformer par étapes un ensemble de vertex, définis dans un espace en trois
dimensions, en un tableau de pixels en deux dimensions.

\begin{figure}[ht]
  \centering
  \includegraphics[width=\textwidth]{figures/openGLPipeline}
  \caption{Vue de haut niveau de la chaîne de traitement de \opengl. Les
    shaders en gras peuvent être définis par l'utilisateur, les shaders en
    pointillé sont optionnels. }
  \label{fig:openGLPipeline}
\end{figure}

Ces étapes peuvent peuvent pour certaines être programmées à l'aide d'un langage
appelé \emph{GLSL} (OpenGL Shading Language), une variante de C, et sont
appelées des \emph{shaders} (voir \textsc{Fig.}~\ref{fig:openGLPipeline}). Voici une
description succinte des shaders programmable utilisés dans \otb, avec leurs
entrées, leurs sorties et une description de leur rôle:
\begin{description}
\item{\bfseries Vertex Shader}
  \begin{itemize}
  \item Entrée: un vertex;
  \item Sortie: un vertex;
  \item Description: ce shader modifie un vertex, par exemple ses attributs de
    position \ç{vec4 gl_Position} ou de taille de point \ç{float gl_PointSize}.
  \end{itemize}
\item{\bfseries Geometry Shader}
  \begin{itemize}
  \item Entrée: un vertex, une ligne, des triangles;
  \item Sortie: zéro, un ou plusieurs vertex, ligne ou triangles;
  \item Description: ce shader permet d'effectuer des calculs géométriques sur
    les primitives qui lui sont données, soit pour modifier leur attributs, soit
    pour générer de nouveaux points, de nouvelles lignes ou de nouveaux
    triangles. Par exemple, dans \otb, ce shader est utilisé pour construire les
    vertex nécessaire à la représentation d'une bactérie, à partir de trois
    vertex donnés en entrée.
  \end{itemize}
\item{\bfseries Fragment Shader}
  \begin{itemize}
  \item Entrée: un fragment;
  \item Sortie: une couleur;
  \item Description: ce shader permet d'attribuer une couleur à un pixel sur
    l'écran à partir des informations calculées par le shader de Rasterization.
    C'est ici par exemple que la représentation de chaque bactérie obtient ses
    couleurs.
  \end{itemize}
\end{description}

\subsubsection{Architecture de \otb}
% Extracted with:
% $ $grep ^module * 2> /dev/null | sed -e 's/^.*\.ml.\?://' | sort -u | copy
L'architecture de \tb{} est organisée en modules \ocaml. L'usage de modules a
été une méthode intuitive pour rapidement prototyper le fonctionnement du
simulateur. Les modules de \tb{} peuvent être découpés en deux catégories: les
modules haut niveau, utilisés dans la programmation de la logique du simulateur
et les modules bas niveau utilisés pour fournir des outils concrets aux modules
de haut niveau, comme par exemple communiquer avec les périphériques
d'entrée/sortie via des binding vers des bibliothèques externes. Les modules
sont écrits avec \ocaml, mais l'accès à \opengl et \opencl ne peut se faire
que via le langage C, ainsi les modules de bas niveau \texttt{Viewer} et
\texttt{OpenCL} sont en partie dédiés à la traduction entre l'interface C et
\ocaml.

\begin{figure}[ht]
  \centering
  \includegraphics[width=\textwidth]{figures/otbModules}
  \caption{Organisation des dépendances entre les modules de \tb{}. Les modules
    de bas niveau sont organisés sur le rang inférieur, les modules de haut
    niveau sont sur le rang supérieur.}
  \label{fig:tbModules}
\end{figure}

\paragraph{Entrée} Le point d'entrée par défaut du programme, à savoir la
fonction d'entrée dans le programme, se trouve dans le module \texttt{Main}. Ce
dernier fait appel aux autres modules pour construire une simulation et
l'exécuter.

\paragraph{Modules de haut niveau}
Les modules de haut niveau regroupent les définitions des fonctions dédié au
travail des modèles: d'une bactérie et d'une population de bactéries, d'un
environnement de morphogènes, d'une zone délimitant les bords physiques de
l'environnement de simulation et enfin de l'assemblage des trois modèles
précédents. Comme indiqué par les flèches de la \textsc{Fig.}~\ref{fig:tbModules}, un
module de haut dépend d'un ou plusieurs modules bas niveau. Le module
\texttt{Bacterium} décrit le modèle d'une bactérie ainsi que l'automate
cellulaire support d'une population de bactéries. Sa fonction la plus notable,
\texttt{Bacterium.run}, prépare et lance le calcul d'un pas de simulation d'une
population de bactéries sur le GPU. Le module \texttt{Morphogen}, assez
symétrique au module \texttt{Bacterium}, décrit l'automate cellulaire modélisant
la réaction-diffusion des morphogènes. Sa fonction principale,
\texttt{Morphogen.run} initialise et lance le calcul d'un pas de simulation
d'une grille de morphogènes sur le GPU. Ces deux modules font appel à quatre
autres modules de bas niveau qui sont décrits dans le paragraphe suivant.
\texttt{Zone} est un module de haut niveau modélisant une boîte de Petri de
taille dynamique en fonction de la répartition des individus qu'elle contient.
Ses bords sont représentées par un cadre jaune dans la
\textsc{Fig.}~\ref{fig:otb2-display}. Finalement, le module \texttt{Coupling} décrit le
modèle complet obtenu par l'interaction entre les bactéries et leur
environnement. C'est ce modèle qui permet aux bactéries de déposer des
morphogènes dans leur environnement et de consommer les morphogènes qui y sont
présent.

\begin{figure}[h]
  \centering
  \includegraphics[width=\textwidth]{otb-multi-view}%
  \caption[Présentation des zones dynamiques]{Trois vues d'une même simulation
    avec \num{1174} bactéries présentant, à gauche, une vue globale montrant
    les zones dynamiques (délimitées en rouge, la zone intérieure est celle
    des bactérie, la zone extérieure est celle des morphogènes). À droite sont
    présentés des vues agrandies de la simulation mettant en évidence, en
    haut, le comportement physique des bactéries et, en bas, la diffusion des
    morphogènes.}
  \label{fig:tbModules}
\end{figure}

\paragraph{Modules de bas niveau}
Le module \texttt{Packfile} regroupe les fonctions de sérialisation nécessaire à
la sauvegarde dans un fichier de l'état d'une simulation. Le module
\texttt{BoundingBox} contient les types et les fonctions permettant de définir
et d'effectuer des calculs courants sur des boîtes encadrantes, une
représentation rectangulaire d'une zone géométrique.

Les modules \texttt{OpenCL} et \texttt{Viewer} sont des modules permettant
d'accéder à des bibliothèques externes. En effet, \tb{} communique avec la
carte graphique \emph{via} le langage et framework \opencl pour les calculs
généraux. Dans un premier temps, nous avons fait appel à un module \ocaml,
appelé SPOC~\cite{bourgoin_spoc:_2012}, fournissant un binding \opencl et une
gestion automatique de la mémoire de la carte graphique. Au delà d'une simple
interface entre \ocaml et \opencl, SPOC abstrait la programmation GPGPU et donne
un accès transparent à d'autres framework. Cependant, SPOC posant des problèmes
de stabilité et possédant bien plus de fonctionnalités que nous avions besoin,
nous avons finalement écrit notre propre interface bas niveau pour \opencl, car
nous souhaitons manipuler explicitement le transfert entre l'hôte et la carte.
Cette interface a pour intérêt d'être simple, stable et de donner une grande
latitude au développeur \ocaml quant à la gestion de la mémoire entre l'hôte et
la carte graphique. Ce binding est contenu dans le module \texttt{OpenCL}. \tb{}
communique aussi avec la carte graphique \emph{via} \opengl. De la même manière,
il existe plusieurs bindings \opengl pour \ocaml. Toutefois, nous n'utilisons
qu'un faible sous-ensemble des fonctions de l'\api{} d'\opengl et nous avons
écrit notre propre interface bas niveau, minimale et spécifique à l'utilisation
de \tb{}. \texttt{Viewer} contient le binding vers \opengl, glew et glfw, trois
bibliothèques utilisées dans la création de l'interface graphique du simulateur.

\subsection{Le langage \sbgp{}}
Pour répondre au point «Entrée» du cahier des charges, nous avons développé
notre propre langage afin de faciliter le paramétrage d'une
simulation \cite{pascalie_developmental_2016} avec \tb{}. \sbgp{} (Synthetic
Biology Genetic Programming) est un langage déclaratif décrivant les options de
configuration d'une simulation \otb. Un fichier de configuration de la
simulation écrit avec \sbgp{} donne une représentation abstraite du génome de la
bactérie sous la forme d'une machine à états finis et inclus des informations
supplémentaires pour paramétrer la simulation, comme les couleurs de rendu.

\paragraph{Syntaxe}
Un fichier \sbgp{} est formaté en JSON~\cite{bray_javascript_2014}. Un génome
\sbgp{} est un objet JSON contenant les blocs suivants:
\begin{itemize}
  \item \texttt{signals}: la liste des signaux en usage dans la simulation,
    chaque signal est définit par un triplet: nom du signal, coefficient de
    diffusion et coefficient de dégradation;
  \item \texttt{reactions}: la liste des réactions, si elles existent. Chaque
    réaction est définit par un triplet: $(r,p,\tau)$ où $r$ est une liste des
    réactants, $p$ une liste des produits et $\tau$ le taux de réaction;
  \item \texttt{type}: la liste des types identifiés par une étiquette unique;
  \item \texttt{behavior}: le comportement, un objet donnant pour chaque type la
    liste des actions à exécuter à chaque pas de temps;
  \item \texttt{transition}: une matrice de taille
    $card(\text{\texttt{type}}) \times card(\text{\texttt{type}})$
  qui représente la table de transition d'état entre chaque types. La matrice
  donne l'identifiant de la transition si il existe.
  \item \texttt{cond\_transition}: la transition conditionnelle identifiée par
    son étiquette; chaque étiquette dans la table de transition d'état doit
    être défini dans ce bloc
\end{itemize}

Le type de chaque bloc est le suivant:
\begin{SBGPcode}
{
  "signals"         : [[Sid_i, Deg_rate_i, Diff_rate_i], [...]],
  "reactions"       : [[[Pid_i, Pid_j], [Rid_i], Rate_i], [...]],
  "type"            : [T0, ..., Tn],
  "behavior"        : {T0 : [{InstructionID : [Pid_1, ..., Pid_n]}]}
  "transition"      : [type x type] of Condition_name,
  "cond_transition" : [{Condition_name_i : ConditionID}, ..., {...}]
}
\end{SBGPcode}

\paragraph{Déclaration des signaux}
Les signaux sont déclarés dans une liste. Chaque signal est un triplet contenant
le nom du signal sous forme de chaîne de caractères, et les coefficients de
diffusion et de dégradation sous forme de nombres réels.
\begin{SBGPcode}
  "signals" : [
    ["alpha", 0.05, 0.01],
    ["beta", 1.5, 1]
  ]
\end{SBGPcode}

\paragraph{Déclaration des réactions}
Les réactions sont déclarés dans une liste. Chaque réaction consiste en un
triplet: la liste des réactants, la liste des produits et la constante cinétique
de la réaction. Par exemple, le bloc \sbgp{} suivant est utilisé pour décrire la
réaction $2\alpha + \beta \xrightarrow{0.15} \gamma$:
\begin{SBGPcode}
  "reactions" : [
    [
      ["alpha", "alpha", "beta"],
      ["gamma"],
      0.15
    ]
  ]
\end{SBGPcode}

\paragraph{Déclaration des types}
Les types de bactéries légaux sont déclarés dans une liste d'identifiants. Le
bloc suivant donne la liste des types impliqués dans le génome:
\begin{SBGPcode}
  "type" : [
    "QUIESCENT",
    "LEADER",
    "FLESH",
    "SKIN",
    "DEAD"
  ]
\end{SBGPcode}

\paragraph{Déclaration du comportement}
Les comportements spécifiques aux types sont déclarés dans une liste
associative; une liste d'instructions est donnée pour chaque étiquette
déclaré dans le bloc type. Par exemple, le code suivant donne le comportement
d'une bactérie de type \texttt{LEADER}:
\begin{SBGPcode}
  "behavior" : {
    "LEADER" : [ { "EmitSignal" : ["alpha", "150"] },
                 { "StopDivide" : [] },
                 { "Accumulate" : ["lambda","1"] }
    ]
  }
\end{SBGPcode}
Une instruction est décrite conformément à la structure suivante:
\begin{SBGPcode}
{ "InstructionID" : ["param1", ...  ,"paramN"] }
\end{SBGPcode}
La valeur du champ \texttt{InstructionID} doit être pris parmi l'ensemble des
instructions \sbgp{} suivantes:
\begin{itemize}
\item \texttt{EmitSignal(s,q)}: dépose une quantité \ç|q| du morphogène \ç|s| à
  la position courante de la bactérie dans l'environnement;
\item \texttt{AbsorbSignal(s,q)}: consomme une quantité \ç|q| du morphogène
  \ç|s| à la position courante de la bactérie dans l'environnement;
\item \texttt{Differentiate(A)}: force un changement de type de la bactérie vers
  le type \ç|A| (il doit exister dans la liste des types);
\item \texttt{Die()}: détruit la bactérie. Elle est en pratique marquée comme
  inactive et réinitialisée;
\item \texttt{Freeze()}: place le marqueur \ç|immobile| sur une bactérie. Dans cet
  état une bactérie est immobile et ne croît pas;
\item \texttt{Unfreeze()}: retire le marqueur \ç|immobile|;
\item \texttt{Divide(r)}: définit le taux de croissance de la bactérie à la
  valeur particulière \ç|r|;
\item \texttt{StopDivide()}: définit le taux de croissance de la bactérie comme
  nul, l'empêchant ainsi de se diviser;
\item \texttt{Accumulate(p,q)}: incrémente le compteur de morphogènes \ç|p| avec
  la valeur \ç|q|;
\item \texttt{Deplete(p,q)}: décrémente le compteur de morphogènes \ç|p| avec la
  valeur \ç|q|.
\end{itemize}

\paragraph{Déclaration des transitions}
Le bloc transition est une matrice encodant le graphe de différentiation du
génome. La taille de la matrice dépend du nombre de types. Lorsque qu'une
transition existe entre deux types, une étiquette est définie dans la case de la
matrice correspondante, sinon, le mot-clé \texttt{NA} est spécifié. Le code
\sbgp{} suivant présente le graphe de différenciation d'un cœur homéostatique à
deux couches: la population croît autour d'une bactérie mère en formant deux
anneaux concentriques constitués d'une quantité invariante de bactéries.
\begin{SBGPcode}
"transition" : [
  ["NA","CL","C1" ,"C2" ,"NA" ,"NA" ],
  ["NA","NA","NA" ,"NA" ,"NA" ,"NA" ],
  ["NA","NA","NA" ,"C12","NA" ,"NA" ],
  ["NA","NA","C21","NA" ,"C23","NA" ],
  ["NA","NA","NA" ,"NA" ,"NA" ,"C3D"],
  ["NA","NA","NA" ,"NA" ,"NA" ,"NA" ]
],
\end{SBGPcode}

\paragraph{Déclaration des conditions de transition}
Une condition est une expression booléenne renvoyant \ç|true| ou \ç|false|. La
valeur du champ \texttt{ConditionID} doit être prise parmi l'ensemble des
instructions \sbgp{} suivantes:
\begin{itemize}
\item \texttt{AND(Cid1,Cid2)}, \texttt{OR(Cid1,Cid2)}: \ç{AND} et \ç{OR} sont les
  opérations classiques de la logique booléenne;
\item \texttt{EqThreshold(s,q)}: ce prédicat n'est vérifié que si la
  concentration du morphogène \ç|s| est égale à \ç|q|;
\item \texttt{LessThreshold(s,q)}: ce prédicat n'est vérifié que si la
  concentration du morphogène \ç|s| est inférieure à \ç|q|;
\item \texttt{GreaterThreshold(s,q)}: ce prédicat n'est vérifié que si la
  concentration du morphogène \ç|s| est supérieure à \ç|q|;
\item \texttt{LessEqThreshold(s,q)}: ce prédicat n'est vérifié que si la
  concentration du morphogène \ç|s| est inférieure ou égale à \ç|q|;
\item \texttt{GreaterEqThreshold(s,q)}: ce prédicat n'est vérifié que si la
  concentration du morphogène \ç|s| est supérieure ou égale à \ç|q|;
\item \texttt{BetweenThreshold(s,q1,q2)}: ce prédicat n'est vérifié que si la
  concentration du morphogène \ç|s| est comprise strictement entre \ç|q1| et
  \ç|q2|;
\item \texttt{InternalProtEqCond(p,q)}: ce prédicat n'est vérifié que si le
  compteur de morphogènes \ç|p| est égal à la valeur \ç|q|;
\item \texttt{InternalProtLessCond(p,q)}: ce prédicat n'est vérifié que si le
  compteur de morphogènes \ç|p| est inférieur strictement à la valeur \ç|q|;
\item \texttt{InternalProtGreatCond(p,q)}: ce prédicat n'est vérifié que si le
  compteur de morphogènes \ç|p| est supérieur strictement à la valeur \ç|q|;
\item \texttt{InternalProtLessEqCond(p,q)}: ce prédicat n'est vérifié que si le
  compteur de morphogènes \ç|p| est inférieur ou égal à la valeur \ç|q|;
\item \texttt{InternalProtGreatEqCond(p,q)}: ce prédicat n'est vérifié que si le
  compteur de morphogènes \ç|p| est supérieur ou égal à la valeur \ç|q|;
\item \texttt{MidPlane(a,b,t)}: ce prédicat n'est vérifié que si la valeur de la
  concentration du morphogène \ç|a| moins la concentration du morphogène \ç|b|
  est comprise entre \ç|-t| et \ç|t|.
\end{itemize}



\section{Propagation Parallèle à la Margolus}\label{sec:ppm}
Nous avons présenté l'architecture de haut niveau de \otb, ainsi que son langage
de configuration, cette section présente l'algorithme de Propagation Parallèle à
la Margolus (\ppm). Cette algorithme forme le socle théorique que nous
implémentons dans les kernels \opencl. Nous commençons par un rappel sur les
automates cellulaires, les automates à blocs, puis nous présentons l'algorithme
de bounding-box dynamique et enfin \ppm.

\subsection{Automates cellulaires et voisinage de Margolus}
% Voir peut-être aussi CML (https://en.wikipedia.org/wiki/Coupled_map_lattice)
Notre travail se situe dans le contexte des automates cellulaires. Un automate
cellulaire est un modèle de calcul discret, introduit par John Von Neumann et
Stanislaw Ulam dans les années 1960~\cite{neumann_theory_1966}. C'est un modèle
de calcul Turing complet qui a d'abord été objet de recherche académique sur les
processus d'auto réplication, ensuite popularisé par le «Game of Life» de John
Conway\footnote{lors de sa publication dans les pages du Scientific American en
  octobre 1970} où il est devenu un modèle pour les processus biologiques. Les
automates cellulaires sont un outil des paradigmes de programmation non
conventionnels, d'un côté permettant de formaliser les processus biologiques
observés, et de l'autre de les modéliser et de les simuler.

\begin{mpo-definition}[Configuration]
  Soit $\oQ$ un ensemble d'états, $G$ un groupe, $c_\oQ : G \rightarrow \oQ$ est
  une \emph{configuration} et $C_\oQ := \{ c \,|\, c:G \rightarrow \oQ \}$ est
  \emph{l'ensemble de toutes les configurations} $c_\oQ$. Quand le contexte lève
  tout ambigüité, on pourra omettre l'indice $\oQ$.
\end{mpo-definition}


\begin{mpo-definition}[Automate cellulaire]\label{def:ca}
  Un automate cellulaire est un tuple $\oA = (M, \oQ, N, f)$ dont:
  \begin{itemize}
  \item $(M,+)$ est un groupe abélien libre appelé \emph{maillage des cellules};
  \item $\oQ$ est \emph{l'ensemble des états}, aussi parfois appelé alphabet,
  \item $N \in M^k$, $k \in \mathbb{N}$ est \emph{l'indice de
      voisinage}, dont les éléments sont uniques deux à deux et $N =
    (n_1,n_2,\ldots,n_k)$;
  \item $f: \oQ^k \rightarrow \oQ$ est la \emph{fonction de transition locale}
  \end{itemize}
  $f$ induit une fonction de transition globale $F$, $\forall c \in C_\oQ$ et $m
  \in M$:
  \[ F(c)(m) := f(c(m+n_1),c(m+n_2),\ldots,c(m+n_k)) \]
\end{mpo-definition}

La mise à jour d'une configuration se fait de manière locale, par l'application
de $f$ sur chaque cellule en fonction du voisinage, et \emph{synchrone} avec
des pas de temps discrets. Dans le cas d'un voisinage classique ou de Von
Neumann, à gauche dans la \textsc{Fig.}~\ref{fig:voisinage-2D}, les voisins directs
sont pris en compte pour le calcul du nouvel état de chaque cellule. Du point de
vue de \mgs{}, un automate cellulaire est constitué d'un ensemble $C$ de
collections topologiques et d'une transformation $T$ jouant le rôle de la
fonction de transition globale de l'automate. Quelques restrictions s'appliquent
sur chaque $c_i \in C$ et $T$: $c_i$ doit être une collection newtonienne de
type \gbf{} et $T$ doit utiliser la stratégie de réécriture
\texttt{maximal-parallel}, afin que la mise à jour soit synchrone.

\begin{figure}[ht]
  \centering
  \includegraphics[page=1,width=.3\textwidth]{figures/ca-neighborhood-2D}
  \includegraphics[page=2,width=.3\textwidth]{figures/ca-neighborhood-2D}
  \includegraphics[page=3,width=.3\textwidth]{figures/ca-neighborhood-2D}
  \caption{Trois exemples de voisinages sur un grille carrée en 2D, à gauche un
    voisinage de Von Neumann, au centre un voisinage hexagonal et à droite un
    voisinage de Moore.}
  \label{fig:voisinage-2D}
\end{figure}

En tant que modèle de calcul, un automate cellulaire fonctionne de manière
fortement parallèle: la même fonction est appliquée sur chaque cellule de la
grille, ou autrement dit le même algorithme est appliqué à de multiples données
ce qui est la définition d'une architecture SIMD. On pourrait donc tirer parti
d'une architecture fortement parallèle pour implémenter un automate cellulaire.
Néanmoins, la dépendance entre les données par la relation de voisinage entre
les cellules pose encore problème. Une solution classique revient à utiliser une
mémoire intermédiaire (un buffer) depuis lequel on va lire les données et écrire
les résultat dans un second buffer, cassant ainsi la dépendance en écriture sur
les cellules concernées. Cependant, il reste une \emph{dépendance en lecture}
sur les cellules voisines causant des collisions de lecture et pouvant encore
porter atteinte à l'efficacité d'une implémentation parallèle. Pour résoudre ces
problèmes, nous nous sommes tournés vers le voisinage de Margolus et de manière
plus générale, le modèle de calcul des automates à blocs.

S'intéressant au développement d'ordinateurs simulant des systèmes physiques,
Toffoli et Margolus introduisent
dans~\cite{toffoli_cellular_1987,margolus_physics-like_1984} les automates
cellulaire à blocs (\emph{block cellular automata}), ou automate à partitions
(\emph{partitioning automata}), pour modéliser le comportement dynamique de
l'écoulement d'un fluide en respectant le principe de \emph{réversibilité} des
lois physiques. En thermodynamique, durant un processus réversible l'entropie
n'augmente pas et le système est en état d'équilibre thermodynamique. On appelle
\emph{automate cellulaire réversible} les automates cellulaire dont on retrouve
la configuration d'origine en «inversant» la flèche du temps. Un automate
cellulaire dont chaque cellule est mise à jour en inversant sa valeur précédente
et en ignorant son voisinage est trivialement réversible. Cependant, au delà
d'une dimension, la réversibilité d'un automate cellulaire est
indécidable~\cite{kari_reversibility_1994}. L'automate gaz sur réseau (lattice
gas automaton)~\cite{hardy_time_1973}, la règle
«Tron»~\cite{toffoli_cellular_1987}, l'ordinateur à boules de billard
(billiard-ball computer)~\cite{fredkin_conservative_1982} sont des exemples
d'automates cellulaires réversibles. Toffoli et Margolus montrent également que
la classe des automates à partitions est équivalente à celles des automates
cellulaires, et que ces deux modèles peuvent se simuler l'un l'autre.

\begin{mpo-definition}[Automate cellulaire à blocs] Un automate cellulaire à
  blocs, aussi appelé automate à partitions, est un tuple
  $B=(M,\oQ,p,\text{\em shift},f)$ avec:
  \begin{itemize}
    % Lattice (order) != Lattice (groups) | Treilli != Réseau
  \item $(M,+)$ le maillage des cellules. C'est le groupe abélien \emph{libre}
    engendré par l'ensemble fini $A = \{g_1,...,g_n\}$;
  \item $\oQ$ un ensemble d'états;
  \item $p : M \rightarrow M$ la fonction de pavage qui pour $\{c_1,...,c_n\}
    \subset \mathbb{Z}^n$ et $\{k_1,...,k_n\} \subset \mathbb{N}_{+}^n$ associe à
    chaque point de $M$ un \emph{bloc}:
    \[
      p(\sum\limits_{i=1}^{n}{c_i\ g_i})
        = \sum\limits_{i=1}^{n}{\left\lfloor \frac{c_i}{k_i} \right\rfloor g_i}
    \]
    La fonction inverse $p^{-1}(m) = \{ m \in M\, |\, p(m) \in M \}$ permet
    d'obtenir l'ensemble des éléments du bloc indicé par $m$.
  \item $\text{\em shift}$ est la liste d'éléments $(m_1,...,m_n)$ de $M$ .
  \item $f : Q^k \rightarrow Q^k$ est la fonction de transition locale, avec $k
    = \prod_{i=1}^{n}{k_i}$.
  \end{itemize}
  La fonction de transition locale $f$ induit une fonction de transition
  globale, dont la définition est décrite par l'algorithme
  \begin{algorithm}
    \ForAll {$s \in \text{shift}$}
    { \ForAll {$m \in M$} {
        $M \gets f (c_1 (T^s \beta(m)))$ \;
      }
    }
  \end{algorithm}
  en posant:
  \begin{itemize}
  \item $T^s$ l'opérateur de décalage sur les groupes Abéliens, $T^s f(g) =
    f(g+s)$ en prenant un groupe Abélien $G$, $g \in G$ et $f : G \rightarrow
    G$;
  \item $\beta(m) = (p^{-1} \circ p) (m)$ l'ensemble des éléments appartenant au
    bloc de la cellule $m$ et
  \item $c_1(B) : M \rightarrow Q^k $ la configuration d'un ensemble de cellules
    \emph{indexées par leur élément correspondant} dans le groupe. Nous
    n'expliciterons pas la version formelle ici pour ne pas alourdir la
    définition inutilement.
  \end{itemize}
\end{mpo-definition}

Nous avons choisi d'utiliser la structure de groupe abélien libre pour
représenter un maillage de cellule. C'est également le point de vue utilisé dans
la définition des \gbf{} à la différence que nous nous restreignons aux groupes
libres, afin de pouvoir définir un pavage naturel. De ce point du vue, il
semblerait qu'un automate cellulaire à bloc ne puisse pas avoir de maillage
hexagonal, avec la présentation de groupe $H = \langle a,b,c \:|\: a + b = c
\rangle$, car ce n'est pas un groupe libre. Néanmoins, nous pourrons considérer
qu'un maillage hexagonal est un maillage carré, de présentation $H = \langle a,b
\rangle$ sur lequel le voisinage entre cellules est spécifique comme présenté
sur la \textsc{Fig.}~\ref{fig:voisinage-2D} au centre. Par conséquent, c'est la
fonction de transition locale qui déterminera quel voisinage est utilisé.

La fonction $\lfloor x/k \rfloor$ se comporte « naturellement » autour de $0$:
chaque image par cette fonction possède exactement $k$ préimages. Cette
caractéristique permet, par exemple, de partitionner $\mathbb{Z}$ en
sous-ensembles de cardinalité $k$, d'où son usage dans la définition de $p$.

Pour la définition de la fonction de pavage, nous avons choisi une approche
similaire à la définition d'une partition de l'espace passant par les
$\mathbb{Z}$-modules comme introduit dans~\cite{michel_representations_1996} à
la section VI.1.1. Notre approche correspond mieux à l'idée intuitive d'une
sous division de l'espace: à chaque cellule du maillage on fait correspondre une
cellule du même maillage, ce qui expose immédiatement la structure et la
régularité du pavage.

Le voisinage de Margolus est la règle de partitionnement du maillage d'un
automate à blocs en une grille composée de blocs de $2$ cellules (ou en carrés
de $2\times 2$ cellules en deux dimensions, en cubes de $2\times 2\times 2$
cellules en trois dimension, …) avec un décalage d'une «case» dans toutes les
directions. C'est donc un cas particulier d'automate cellulaire à blocs, avec
les restrictions suivantes (en reprenant le contexte de la définition des
automates à blocs):
\begin{itemize}
\item $p : M \rightarrow M$ la fonction de pavage qui pour $\{c_1,...,c_n\}
  \subset \mathbb{Z}^n$ associe à chaque point de $M$ un bloc de côté 2:
    \[
      p(\sum\limits_{i=1}^{n}{c_i\ g_i})
        = \sum\limits_{i=1}^{n}{\left\lfloor \frac{c_i}{2} \right\rfloor g_i}
    \] et
\item shift est l'élément $\sum\limits_{i=1}^{n}{1\ g_i} \in M$.
\end{itemize}

\begin{figure}[ht]
  \centering
  \includegraphics[page=1,width=.4\textwidth]{figures/ca-neighborhood}
  \hspace{1cm}
  \includegraphics[page=2,width=.4\textwidth]{figures/ca-neighborhood}
  \caption{Voisinage classique (à gauche) d'un automate cellulaire à une
    dimension et voisinage de Margolus (à droite). Pour le voisinage classique,
    chaque ligne correspond à un pas dans la simulation de l'automate, de haut
    en bas. Pour le voisinage de Margolus, un seul pas de simulation est
    présenté, la ligne du milieu est présentée afin de montrer la correspondance
    entre les deux automates.}
  \label{fig:ac-nh}
\end{figure}


\subsection{Bounding box et activité}\label{sec:boundingbox}
Dans OTB, l'environnement dans lequel évoluent les bactéries n'a pas \emph{a
  priori} de taille fixe et subit une expansion ou une contraction en fonction
de l'espace occupé par les bactéries, ou par leur support. Nous appelons
\emph{bounding box} la boîte englobante dont l'intérieur forme l'environnement
dans lequel évoluent les bactéries, ou leur support. Cette bounding box a une
taille rectangulaire libre. La bounding box est construite \emph{a posteriori}:
elle est déterminée par la taille de la grille d'automate qu'elle représente. De
cette manière, il est possible d'accéder rapidement à la taille d'une population
de bactéries et d'effectuer des calculs de collision entre bounding box et des
opérations de séparation et de fusion. Ces opérations, bien que planifiées, ne
sont que partiellement implémentées à l'heure de l'écriture de ce mémoire.

\begin{figure}[ht]
  \centering
  \includegraphics[page=1,width=.8\textwidth]{figures/bb-activity}
  \caption{Grille d'un automate cellulaire, ses indicateurs vertical et
    horizontal et sa Bounding Box (en gras)}
  \label{fig:bb-activity}
\end{figure}

Nous utilisons une mesure d'activité spatiale pour déterminer s'il est
nécessaire d'agrandir ou de diminuer la taille de l'environnement. Chaque
cellule de l'automate contenant une bactérie est considérée comme active. Si une
cellule active se trouve à une distance inférieure à $d_c$ d'un des bords, alors
l'environnement est agrandi. À l'inverse, si la distance la cellule active la
plus proche d'un bord à ce bord est supérieure à $d_f$ alors, l'environnement
est rétréci.

En pratique, une bounding box est définie par un triplet $(p,w,h)$ où $p=(x,y)$
est un point dans le plan représentant le coin inférieur gauche de la boîte, $w$
est sa largeur et $h$ est sa hauteur (Fig.~\ref{fig:bb-activity}). Nous mesurons
l'occupation pour chaque ligne et chaque colonne avec deux tableaux appelés
\emph{indicateur}: un pour l'occupation horizontale, l'autre pour l'occupation
verticale. Ces deux tableaux sont suffisant pour déterminer s'il est nécessaire
d'agrandir ou de rétrécir l'environnement. L'algorithme suivant est utilisé par
\otb pour déterminer une nouvelle bounding box optimale, $w,h$ sont la largeur
et la hauteur respectives de la bounding box de départ, $c_w,c_h$ sont les
tableaux indicateurs sur la largeur et la hauteur respectivement.

La première partie de cet algorithme détermine la position de coupe par rapport à
chacun des quatre bords nord, sud, est et ouest de la bounding box. Il suffit
pour cela de déterminer la position du premier élément non nul pour chaque
extrémité de chaque indicateur, ce qui nous fournit le quadruplet $(\text{posN},
\text{posS}, \text{posE}, \text{posW})$ de positions de coupe au nord, au sud à
l'est et à l'ouest respectivement.

La seconde partie de cet algorithme détermine la nouvelle dimension de la
bounding box, en prenant en compte les deux marges notée $d_c$ et $d_f$ et en
fonction des positions de coupe déterminées auparavant. La fonction
\texttt{CropBoundingBox()} est chargée de la mise en œuvre du découpage de
l'ancienne bounding box et en retourne une nouvelle en fonction du quadruplet
$(\text{n}, \text{s}, \text{e}, \text{w})$ indiquant les corrections de marge à
chacune des extrémité. Par exemple, $\text{n}=2$ agrandi de 2 unités la taille
de la bounding box au nord, alors que $\text{e}=-4$ rétrécie de 4 unités la
taille de la bounding box à l'est. Cette fonction est implémenté dans le module
de bas niveau \ç|BoundingBox|.

\begin{algorithm}[H]
  \SetKwInOut{Input}{input}\SetKwInOut{Output}{output}
  \SetKwData{PosN}{posN}\SetKwData{PosS}{posS}
  \SetKwData{PosE}{posE}\SetKwData{PosW}{posW}
  \SetKwData{OldBB}{oldBB}
  \SetKwFunction{CBB}{cropBoundingBox}

  \Input{Bounding box actuelle \OldBB et \PosN, \PosS, \PosE, \PosW}%
  \Output{Une nouvelle bounding box ajustée}
  \eIf(\tcp*[f]{Marge nord}){$\PosN < d_c$}{
    $n \gets d_c - \PosN$\;
  }{
    \If{$\PosN > d_f$} {
      $n \gets d_c - \PosN$\;}
    \Else{
      $n \gets 0$\;
    }
  }

  \tcp{Exactement la même opération pour sud, est et ouest}

  \Return{\CBB{\OldBB, $n$, $s$, $e$, $w$}}
  \caption{Création d'un nouvelle bounding box aux \emph{bonnes} dimensions}
\end{algorithm}



\subsection{Propagation Parallèle à la Margolus}
L'algorithme \ppm est inspiré du voisinage de Margolus, à droite dans la
\textsc{Fig.}~\ref{fig:ac-nh}, en deux dimensions et l'utilise comme structure de
données. Il effectue la simulation d'un automate cellulaire à blocs sur une
machine parallèle d'architecture SIMD, nous l'avons implémenté sous forme d'un
kernel \opencl. Pour faciliter sa lisibilité, en voici une description un peu
plus haut niveau que notre implémentation.


\paragraph{Cellule}%
Chaque cellule est constituée de quatre emplacements: son propre état et l'état
de trois de ses voisines comme présenté dans la \textsc{Fig.}~\ref{fig:margolus-grid}.
Ces quatre emplacements seront notés dans l'algorithme suivant par abus
de la notation classique de tableau, avec les indices spéciaux 1, 2, 3 et
4. Ainsi, $m[a1]$ désigne le premier emplacement de la cellule $A$. La
\textsc{Fig.}~\ref{fig:margolus-grid} (a) donne une représentation visuelle de la
position des cellules et de leurs emplacements en accord avec l'algorithme.

\begin{figure}[ht]
  \centering
  \includegraphics[page=1,width=.24\textwidth]{figures/margolus}\hfill%
  \includegraphics[page=2,width=.24\textwidth]{figures/margolus}\hfill%
  \includegraphics[page=3,width=.24\textwidth]{figures/margolus}\hfill%
  \includegraphics[page=4,width=.24\textwidth]{figures/margolus}

  \hspace{.11\textwidth}(a)\hfill(b)\hfill(c)\hfill(d)\hspace{.11\textwidth}
  \caption[Quatre configuration successives d'un automate cellulaire en deux
    dimensions.]{\textit{Général.} Quatre configurations successives d'un
    automate cellulaire en deux dimensions dont les cellules sont en gris.%
    \textit{\hspace{1em}(a).} Les coordonnées utilisées dans l'algorithme pour
    une cellule de Margolus.%
    \textit{\hspace{1em}(b).} Les cellules A,B,C et D font partie d'une cellule
    du voisinage de Margolus (en noir, en gras).%
    \textit{\hspace{1em}(c).} Chacune des cellules stocke quatre états: la copie
    de l'état de trois de ses voisines et le sien (en blanc sur fond noir).
    Le positionnement de ces états est présenté par rapport au voisinage de
    Margolus courant.%
    \textit{\hspace{1em}(d).} Après une étape de calcul, les états de A, B, C
    et D sont mis à jour en A', B', C' et D' et le voisinage de Margolus est
    décalé.}
  \label{fig:margolus-grid}
\end{figure}

\paragraph{Algorithme}%
\ppm prend en entrée un tableau de cellules \ç|cellSpace[w][h]| en deux
dimensions de largeur \ç|w| et de hauteur \ç|h| et un indicateur \ç|shift|
déterminant le décalage du voisinage. La fonction d'évolution locale
\ç|evol(c,n1,n2,n3,n4,n5,n6,n7,n8)| prend en argument l'état courant \ç|c| ainsi
que l'état des 8 cellules voisines \ç|n1,...,n8|.

La fonction de mise à jour est appelée en parallèle sur \emph{le même tableau}
de cellules et doit donc être paramétrée pour indiquer sur quelle cellule elle
travaille: c'est le rôle des deux variables spéciales $gc_x$ et $gc_y$ indiquant
la position globale en abscisse et en ordonnée du bloc courant. C'est donc la
coordonnée d'une cellule du voisinage de Margolus qui est fournie. Ces deux
variables sont les paramètres propres à un work-item \opencl, qui permet à
chaque instance du kernel d'avoir un comportement différent.

\begin{algorithm}
  \SetKwInOut{Input}{input}\SetKwInOut{Output}{output}
  \SetKwData{CellSpace}{cellSpace}\SetKwData{Shift}{shift}
  \SetKwData{Slot}{.slot}\SetKwData{Evol}{evol}

  \Input{Une configuration \CellSpace d'un automate cellulaire en deux
    dimension, un indice de décalage \Shift valant 0 ou 1 et les coordonnées
    globales $(gc_x,gc_y)$ du bloc actuel}%
  \Output{Aucune (le changement est effectué par effet de bord)}

  $p_x \gets 2 * gc_x + \texttt{shift}$ \tcp*{Coordonnées de la cellule}
  $p_y \gets 2 * gc_y + \texttt{shift}$\;
  $i_x \gets p_x   + 1$\;
  $i_y \gets p_y   + 1$\;

  Nouveau tableau local de 2×2 cellules: $m[2][2]$\;
  $m[0][1] \gets \texttt{cellSpace}[N*iy+px]$\;
  $m[1][1] \gets \texttt{cellSpace}[N*iy+ix]$\;
  $m[0][0] \gets \texttt{cellSpace}[N*py+px]$\;
  $m[1][0] \gets \texttt{cellSpace}[N*py+ix]$\;

  Nouvelle cellule locale: $n$\;
  $n[1] \gets \Evol(m[a3],$\tcp*[f]{Cellule courante $A$}\\
  $\qquad m[a1],m[a2],m[b1],m[b4],m[c1],m[d2],m[d1],m[a4])$
  \;

  $n[2] \gets \Evol(m[b4],$\tcp*[f]{Cellule courante $B$}\\
  $\qquad m[a2],m[b1],m[b2],m[b3],m[c2],m[c1],m[d2],m[a3])$
  \;

  $n[3] \gets \Evol(m[c1],$\tcp*[f]{Cellule courante $C$}\\
  $\qquad m[a3],m[b4],m[b3],m[c2],m[c3],m[c4],m[d3],m[d2])$
  \;

  $n[4] \gets \Evol(m[d2],$\tcp*[f]{Cellule courante $D$}\\
  $\qquad m[a4],m[a3],m[b4],m[c1],m[c4],m[d3],m[d4],m[d1])$
  \;

  $\CellSpace[N*iy+px] \gets n$\;
  $\CellSpace[N*iy+ix] \gets n$\;
  $\CellSpace[N*py+px] \gets n$\;
  $\CellSpace[N*py+ix] \gets n$\;
  \caption{\ppm}
\end{algorithm}

La cellule localement créée $n$ est au cœur du fonctionnement de \ppm.
Il pourrait paraître en premier lieu étrange qu'elle soit recopiée à
l'identique dans les quatre cellules de \ç|cellSpace| constituant l'unique
cellule de Margolus en cours de mise à jour. Néanmoins, en s'aidant de la
\textsc{Fig.}~\ref{fig:margolus-grid} et en prenant en compte que le prochain
décalage de la grille de Margolus se fait sur une cellule en diagonale, une
symétrie apparaît.

La création d'un tableau local $m$ de $2 \times 2$ est en revanche moins
importante, cependant elle reste utile pour deux raisons: elle améliore sa
lisibilité de l'algorithme et l'utilisation de variables locales en général nous
permet de limiter le nombre d'accès à la mémoire globale en faisant un usage
plus intensif de la mémoire privée (voir le paragraphe Architecture \opencl
page~\pageref{sec:archi-cl}).

L'algorithme \ppm travaille sur un bloc et, dans le voisinage de Margolus,
chaque bloc est indépendant, ainsi il n'y a pas de problème de concurrence ni à
l'accès au données ni à l'écriture. Chaque cellule embarque un quart de son
voisinage et comme un bloc est constitué de quatre cellules, ce bloc contient
toutes les informations nécessaires pour effectuer une mise à jour localement.
Ainsi, alors que les algorithmes usuels de simulation d'automates cellulaires
utilisent deux copies de la configuration actuelle pour éviter les collisions de
lecture, \ppm n'en a pas besoin.







% #IloveFloatBarrier
\FloatBarrier
\section{Implémentation de \tb{}}\label{sec:implementation}
Dans cette section nous décrivons les trois modèles (appelés moteurs) qui
gouvernent le fonctionnement du simulateur: le moteur chimique, le moteur
physique et le moteur de décision.

\subsection{Moteur chimique: réaction et diffusion}
Dans cette section nous présentons la fonction de transition de l'automate
cellulaire en charge de la simulation chimique du milieu dans lequel évolueront
les bactéries. Ce milieu a comme propriété de diffuser plusieurs morphogènes et
de les faire réagir. L'équation de réaction-diffusion est un modèle mathématique
qui permet, entre autres, de modéliser l'évolution de la concentration dans
le temps et dans l'espace d'espèces chimiques diffusant et réagissant entre
elles. Il a notamment été utilisée par A. Turing~\cite{turing_chemical_1952}
a qui nous empruntons le terme de \emph{morphogène} pour décrire une espèce
chimique dans notre automate. Exprimée pour un seul morphogène dans une
dimension d'espace, l'équation de réaction-diffusion est appelé l'équation de
Kolmogorov-Petrovsky-Piskounov (KPP)~\cite{kolmogorov_etude_1937}.

Nous partirons d'une version généralisée de l'équation KPP en deux dimensions
d'espace, avec un nombre $n \in \mathbb{R}_+$ de morphogènes:
\begin{equation}
  \frac{\partial\orr\varphi}{\partial t} =
  \orr{\text{X}}(\orr\varphi) +
  \orr{\text{D}}^{\text{T}} \cdot \left( \frac{\partial^2\orr\varphi}{\partial x^2} +
                                  \frac{\partial^2\orr\varphi}{\partial y^2} \right)
  \label{eq:reaction-diffusion-2D}
\end{equation}
L'équation~\ref{eq:reaction-diffusion-2D} donne la variation des concentrations
de plusieurs espèces chimiques, ou morphogènes, dénotée $\orr\varphi =
(\varphi_1,\ldots,\varphi_n)$, au cours du temps et dans un espace à deux
dimensions où $\orr{\text{D}}^{\text{T}} = (\text{D}_1,\ldots,\text{D}_n)$ sont
les coefficients de diffusion et $\orr{\text{X}}^{\text{T}} =
(\text{X}_2,\ldots,\text{X}_n)$ est un champ de vecteur décrivant la cinétique
locale du système, autrement dit les réactions entre morphogènes.

\paragraph{Discrétisation}
Lors d'une simulation numérique, le temps et l'espace sont discrétisés et il
existe différentes méthodes pour approcher par des méthodes discrètes la
solution de l'équation~\ref{eq:reaction-diffusion-2D}. L'implémentation actuelle
du simulateur utilise la méthode d'Euler, la méthode d'intégration numérique la
plus simple:
\begin{align}\label{eq:num-int}
  \displaystyle
  \varphi^{t + \Delta t}_{k,i,j} &= \Delta t\,\text{X}_1(\varphi^t_{1,i,j}, \ldots, \varphi^t_{k,i,j})
  + \frac{\text{D}_k\Delta t}{(\Delta x)^2}
  (\varphi^t_{k,i-1,j} + \varphi^t_{k,i+1,j} + \varphi^t_{k,i,j-1} + \varphi^t_{k,i,j+1}
  - 4\varphi^t_{k,i,j})
\end{align}
où les indices $(i,j)$ dénotent la position sur une grille carrée en 2D.
Dans~\cite{dilao_validation_1998}, les auteurs fournissent une classe
de méthodes pour l'intégration numérique d'équations de diffusion et de
réaction-diffusion approchant la solution exacte de l'équation de diffusion
continue. En l'absence d'une validation plus poussée de notre outil, nous ne
pouvons pas encore nous prononcer sur l'impact de l'erreur introduite par la
discrétisation sur automate. Cependant, d'un point de vue qualitatif, le moteur
chimique s'est suffisamment bien comporté pour les quelques tests que nous avons
menés. En choisissant un $\Delta t$ suffisamment petit, en l'occurrence $\Delta
t = \num{0.01}$, nous avons pu reproduire des réactions chimiques classiques,
dont une réaction BZ (voir section~\ref{sec:bz}).


\paragraph{Fonction de transition}
\newcommand{\RTable}{\ensuremath{\mathit{RTable}}}
\newcommand{\RRate}{\ensuremath{\mathit{RRate}}}
Chaque cellule de l'automate a pour état la concentration de chacun des $n$
morphogènes de la simulation:
\[
  (conc_1,conc_2,\ldots,conc_n)
\]
Une table des réactions $\RTable[i][j]$ est définie comme un tableau en deux
dimensions avec selon $i$ les $n$ morphogènes et selon $j$ les $r$ réactions.
Une liste des taux de réaction $\RRate[j]$ donne le coefficient de vitesse pour
la réaction $j$. Ces deux valeurs sont définies avant la simulation et sont donc
accessibles à la fonction de transition locale du moteur chimique. Par exemple,
dans le cas d'une simulation à trois morphogènes $a,b,c$ où deux réactions $r,s$
sont possibles
\[
  \RTable = \left[
    \begin{array}{rrr}
      -1 & -2 & 2  \\
      0  & 2  & -1
    \end{array} \right]
  \qquad
  \RRate = (\alpha, \beta)
\]
désigne les réactions
\[
  a + 2b \xrightarrow{\alpha} 2c
  \qquad\text{et}\qquad
  c      \xrightarrow{\beta}  2b
\]

\begin{algorithm}
  \SetKwData{Evol}{evol}
  \SetKw{KwFrom}{from}
  \SetKw{And}{and}
  \SetKwData{NBPC}{$\mathit{nbpc}$}
  \SetKwProg{Fn}{Fonction}{ as}

  \Fn{\Evol(c: cell, $n_1$: cell, $\ldots$, $n_8$: cell): cell}{
    \For(\tcp*[f]{Diffusion de chacun des morphogènes}){$i$ \KwFrom $0$ \KwTo
      $n-1$}{%
      $ c'[i] \gets c[i] + \Delta t * D[i] * \frac{1}{(\Delta x)^2} * ( n_2[i] +
      n_4[i] + n_6[i] + n_8[i] - 4 c[i] )$\;%
      \For(\tcp*[f]{Résultat des réactions entre morphogènes}){$j$ \KwFrom $0$
        \KwTo $r-1$}{ $c'[i] \gets c'[i] + \Delta t * \RTable[i][j] * \RRate[j]$\;
      } } \Return{$c'$}\; }
  \caption{Fonction de transition locale du moteur chimique}
\end{algorithm}

Dans le cas de la diffusion sans réaction, la diffusion sur
automate à bloc à voisinage de Margolus a été traité dans la
littérature~\cite{medvedev_multi-particle_2010} et montre que le processus
de diffusion est exactement simulé par un automate cellulaire booléen
bloc synchrone sur voisinage de Margolus.


\subsection{Moteur physique: collisions}
Dans cette section, nous présentons la construction de la fonction de transition
de l'automate cellulaire en charge de la simulation physique des bactéries. Le
moteur physique représente certains aspects de la mécanique du solide, notamment
les translations sur le plan, la rotation et la collision entre objets. Nous
utilisons un modèle de mécanique newtonien avec des simplifications spécifiques
à l'échelle à laquelle les bactéries évoluent. Le modèle suivant est décrit dans
un plan vectoriel en deux dimensions.

\paragraph{Modèle physique au repos}
Une bactérie est modélisée par une masse $m$, un centre de masse $p =
(p_x,p_y)$, une longueur $l$ et un rayon $r$ donnant son épaisseur (voir
Fig.~\ref{fig:ecoli} droite). Son orientation est représentée de manière
équivalente, soit par un angle $t$, soit par un vecteur directeur $d =
(d_x,d_y)$ en prenant $d_x = \mathit{cos}(t)$ et $d_y = \mathit{sin}(t)$. La
longueur totale d'une bactérie est notée $L$ et sa surface $S$.
Nous ajoutons deux points $f=(f_x,f_y)$ et $b=(b_x,b_y)$ représentant respectivement l'avant et
l'arrière de la bactérie
\begin{align*}
  f_x &= p_x + l \times \mathit{cos}(t) & b_x &= p_x - l \times \mathit{cos}(t)\\
  f_y &= p_y + l \times \mathit{sin}(t) & b_y &= p_y - l \times \mathit{sin}(t)
\end{align*}
%
% Forces
Du point de vue cinétique, une bactérie possède une vitesse de déplacement de
norme $v_b$ colinéaire à $\orr{d}$ et une vitesse de rotation de norme
$\omega_b$ autour de $p$. Elle a également un taux de croissance $g$ faisant
varier sa masse, sa longueur et sa largeur.

La quantité de mouvement, aussi appelé moment linéaire, d'une bactérie est
$\orr{p} = m \orr{v}_b$ et son moment angulaire, aussi appelé moment cinétique,
autour de son centre d'inertie est $\orr{L} = I \orr{\omega}_b$ où le moment
d'inertie $I$ est similaire au moment d'inertie d'une tige de densité de masse
uniforme, tournant autour de son centre de masse, de la même longueur que la
bactérie, à savoir:
\[
  I = \displaystyle\frac{1}{12} \times m(2r+2l)^2 =
  \displaystyle\frac{1}{6}  \times m(r+l)^2
\]

Du point de vue dynamique, nous considérons que les bactéries évoluent dans un
milieu dont la viscosité est très élevée. Au centre de masse de chaque bactérie
s'applique la somme des forces extérieures d'accélération linéaire $F$ et
accélération angulaire $T$ transmises lors de collisions avec d'autres objets.

% Calcul d'une collision (schéma)
\begin{figure}[ht]
  \centering
  \includegraphics[page=2,width=.49\textwidth]{figures/collision}\hfill
  \includegraphics[page=3,width=.49\textwidth]{figures/collision}
  \caption{La détection d'une collision a lieu lorsque l'un des pôles d'une
    bactérie est proche du corps d'une autre bactérie (à gauche). Quelques
    vecteurs utilisés dans le calcul de l'impulsion (à droite)}
  \label{fig:collision}
\end{figure}

% Qu'est-ce-qu'une collision
\paragraph{Collision entre deux bactéries}
En général, déterminer si deux objets sont en collision, c'est à dire déterminer
si l'intersection des deux objets est non vide, n'est pas trivial, et ce pour
deux raisons: les objets ont une forme potentiellement complexe et il peut
y avoir un grand nombre d'objets. Pour les ceux qui ont une forme complexe,
un objet de forme proche et plus simple, comme un disque ou un rectangle,
appelé \emph{bounding box}\footnote{différent de la bounding box que nous avons
présenté avant, celle-ci est autour d'une unique bactérie, notre bounding
box est autour d'une population de bactérie. Leurs fonctions sont toutefois
similaires.} est utilisée pour faciliter le calcul d'intersection. Dans le cas
où un grand nombre d'objet peuvent entrer en collision, l'espace est divisé en
petites sections afin que le test de collision n'ait lieu que pour un faible
nombre d'objets.

% Collision de deux bactéries
Commençons par déterminer la fonction de collision pour deux bactéries. Cette
fonction prend en entrée deux bactéries et retourne vrai s'il y a collision,
avec la distance de collision la plus courte, les vecteurs de collision, les
positions des collisions, ou faux s'il n'y a pas collision. Chaque bactérie est
déterminée par un quadruplet comprenant le vecteur position du centre de masse,
le vecteur directeur (unitaire), la longueur et le rayon respectivement $b_1 =
(\orr{p_1},\orr{d_1},l_1,r_1)$ et $b_2 = (\orr{p_2},\orr{d_2},l_2,r_2)$. À
partir de ces informations, nous pouvons obtenir les coordonnées des extrémités
\begin{align*}
    \orr{f_1} &= \orr{p_1} + l_1 \cdot \orr{d_1} & \orr{f_2} &= \orr{p_2} + l_2 \cdot \orr{d_2}\\
    \orr{b_1} &= \orr{p_1} - l_1 \cdot \orr{d_1} & \orr{b_2} &= \orr{p_2} - l_2 \cdot \orr{d_2}
\end{align*}
Le critère de collision est simple. Grâce à la forme de la bactérie, il nous
faut juste vérifier qu'aucun point du segment $[f_1,b_1]$ n'est à une distance
inférieure ou égale à $d_\mathrm{min} = r_1+r_2$ d'un point du segment
$[f_2,b_2]$. Aussi, nous commencerons par construire les projections
orthogonales restreintes à un segment des pôles $\orr{f}^j_i, \orr{b}^j_i$ d'une
bactérie $i$ sur le corps d'une bactérie $j$.
\begin{align*}
  \orr{f}^2_1 &= \orr{p_2} +
                  \big[ \orr{d_2} + (\orr{f_1} - \orr{p_2}) \big]^{l_2}_{-l_2}
                  \cdot{} \orr{d_2} &
  \orr{f}^1_2 &= \orr{p_1} +
                  \big[ \orr{d_1} + (\orr{f_2} - \orr{p_1}) \big]^{l_1}_{-l_1}
                  \cdot{} \orr{d_1} \\
  \orr{b}^2_1 &= \orr{p_2} +
                  \big[ \orr{d_2} + (\orr{b_1} - \orr{p_2}) \big]^{l_2}_{-l_2}
                  \cdot{} \orr{d_2} &
  \orr{b}^1_2 &= \orr{p_1} +
                  \big[ \orr{d_1} + (\orr{b_2} - \orr{p_1}) \big]^{l_1}_{-l_1}
                  \cdot{} \orr{d_1}
\end{align*}
Restreindre la projection a un segment reflète le fait qu'une bactérie a une
longueur \emph{limitée}. Si on omet cette restriction, on obtient les
projections de la \textsc{Fig.}~\ref{fig:collision}, où la distance la plus courte
n'est pas nécessairement indicatrice d'une collision. Il nous reste à obtenir
les 4 vecteurs de collisions
\begin{align*}
  \orr{n}_{f_{12}} &= \orr{f}^2_1 - \orr{f}_1 & \orr{n}_{f_{21}} &= \orr{f}^1_2 - \orr{f}_2 \\
  \orr{n}_{b_{12}} &= \orr{b}^2_1 - \orr{b}_1 & \orr{n}_{b_{21}} &= \orr{b}^1_2 - \orr{b}_2
\end{align*}
dont les normes nous donnent bien la distance entre les corps des bactéries. Si
$\| \orr{n}_{f_{12}} \| \le d_\mathrm{min}$ ou $\| \orr{n}_{b_{12}} \| \le
d_\mathrm{min}$ ou $\| \orr{n}_{f_{21}} \| \le d_\mathrm{min}$ ou $\|
\orr{n}_{b_{21}} \| \le d_\mathrm{min}$, alors il existe au moins une collision
entre les deux bactéries. Le point d'impact se trouve sur le segment entre
le pôle de la première bactérie et son projeté sur le corps de la seconde.
Supposons que la collision ait lieu lorsque les corps des deux bactéries sont au
contact, alors les points de collisions peuvent être les suivants
\begin{align*}
  \orr{c}_{f_{12}} &= \displaystyle\frac{r_1 \cdot\orr{f}_1 + r_2 \cdot\orr{f}^2_1}{r_1 + r_2} & \orr{c}_{f_{21}} &= \displaystyle\frac{r_2 \cdot\orr{f}_2 + r_1 \cdot\orr{f}^1_2}{r_1 + r_2} \\
  \orr{c}_{b_{12}} &= \displaystyle\frac{r_1 \cdot\orr{b}_1 + r_2 \cdot\orr{b}^2_1}{r_1 + r_2} & \orr{c}_{b_{21}} &= \displaystyle\frac{r_2 \cdot\orr{b}_2 + r_1 \cdot\orr{b}^1_2}{r_1 + r_2}
\end{align*}



\paragraph{Modèle physique post-collision}
Pour calculer le résultat des échanges de forces lors d'un impact entre une
bactérie $i$ et une bactérie $j$, il est nécessaire de déterminer le type
d'impact que nous souhaitons représenter. Dans le cas d'un impact parfaitement
élastique, toute l'énergie transmise à la bactérie est restituée et l'énergie
cinétique du système est conservée. Dans le cas d'un impact inélastique, une
partie de l'énergie est dissipée et seule une partie de l'énergie totale est
restituée. En une dimension, les vitesses post-collision de deux objets $a$ et
$b$ sont données par les deux équations suivantes
\[
  v_a = \displaystyle\frac{Cr\ m_b (u_b - u_a) + m_a u_a + m_b u_b}{m_a + m_b}
  \qquad
  v_b = \displaystyle\frac{Cr\ m_a (u_a - u_b) + m_a u_a + m_b u_b}{m_a + m_b}
\]
où
\begin{itemize}
\item $Cr \in [0,1]$ est le \emph{coefficient de restitution}, $0$ pour une
  collision parfaitement inélastique et $1$ pour une collision parfaitement
  élastique;
\item $u_a, u_b$ sont les vitesses avant collision des objets $a, b$
  respectivement;
\item $m_a, m_b$ sont les masses des objets $a, b$ respectivement.
\end{itemize}

Dans le cas d'une collision en 2D, nous pourrions appliquer ces deux équations
sur chacune des deux composantes des vecteurs vitesse, mais il nous manquerait
les informations sur le moment angulaire de chaque bactérie. Nous utilisons un
modèle un peu différent, ajusté à la collision de corps en 2D dont voici la
description.

\noindent
Les vecteurs vitesse \emph{au point de collision} pré-collision sont
\[
  \orr{v}_{ac1} = \orr{v}_1 + \orr{\omega}_1 \wedge \orr{r}_1
  \qquad
  \orr{v}_{ac2} = \orr{v}_2 + \orr{\omega}_2 \wedge \orr{r}_2
\]
Notre objectif est de trouver les nouvelles valeurs $\orr{v}'_i$ et
$\orr{\omega}'_i$ des deux bactéries après la collision. Écrivons d'abord les
vecteurs vitesse au point de collision post-collision
\[
  \orr{v}_{pc1} = \orr{v}'_1 + \orr{\omega}'_1 \wedge \orr{r}_1
  \qquad
  \orr{v}_{pc2} = \orr{v}'_2 + \orr{\omega}'_2 \wedge \orr{r}_2
\]
Les vitesses relatives au point d'impact $\orr{v}_{ac12}, \orr{v}_{pc12}$
pré-collision et post-collision respectivement sont données par
\[
  \orr{v}_{ac12} = \orr{v}_{ac1} - \orr{v}_{ac2}
  \qquad
  \orr{v}_{pc12} = \orr{v}_{pc1} - \orr{v}_{pc2}
\]
Dans notre modèle, nous faisons l'hypothèse que la vitesse à laquelle les deux
bactéries se séparent est proportionnelle à la vitesse à laquelle les deux
bactéries se rapprochent.
\begin{equation*}
  \orr{v}_{pc12} \cdot{} \orr{n} = -\mathit{Cr}\ \orr{v}_{ac12} \cdot{} \orr{n}
\end{equation*}
Nous retrouvons ici le coefficient de restitution $\mathit{Cr} \in [0,1]$ et
\orr{n} est le \emph{vecteur normal unitaire} ($\|\orr{n}\| = 1$) au point de
collision pointant vers l'extérieur de la seconde bactérie.

La collision en elle-même est modélisée par une \emph{impulsion}. Pour une
bactérie quelconque, rencontrer un obstacle équivaut à se voir appliquer une
force au point de collision, dont la durée d'application est suffisamment courte
pour qu'elle soit considérée \emph{constante}. Cette application s'effectue
\emph{sans frottements} ce qui implique qu'elle soit dirigée selon \orr{n}, la
normale de collision. Une impulsion $\orr{J} = j \cdot \orr{n}$ est
l'intégration de cette force sur $\Delta_t = t' - t$ le temps d'application
\[
  \orr{J} = \int^{t'}_t \orr{F} dt
\]
Comme d'après la seconde loi de Newton, une force est la dérivée temporelle
d'une quantité de mouvement et que l'on suppose $\Delta_t$ suffisamment court
pour que \orr{F} ne varie pas, une impulsion est une variation entre les
quantités de mouvement entre les temps $t$ et $t'$ et a donc la dimension d'une
quantité de mouvement. La collision modifie les vitesses des deux bactéries
selon \orr{n} avec l'intensité $j$. En divisant le tout par la masse, on
retrouve la dimension d'une vitesse. L'impulsion s'applique selon \orr{n} pour
la première bactérie et selon $-\orr{n}$ pour la seconde.
\begin{equation*}
  \orr{v}'_{1} = \orr{v}_{ac1} + j \orr{n} / m_1
  \qquad
  \orr{v}'_{2} = \orr{v}_{ac2} - j \orr{n} / m_2
\end{equation*}
Les vitesses angulaires sont aussi modifiées de manière symétrique en prenant en
compte la distance au centre de rotation des bactéries — si cette distance est
nulle, la vitesse angulaire est nulle — et en divisant cette fois
par le moment d'inertie de chacune des bactéries
\begin{equation*}
  \orr{\omega}'_{1} = \orr{\omega}_{ac1} + (\orr{r}_1 \wedge j \orr{n}) / I_1
  \qquad
  \orr{\omega}'_{2} = \orr{\omega}_{ac2} - (\orr{r}_2 \wedge j \orr{n}) / I_2
\end{equation*}

\noindent
Après résolution des équations précédentes pour $j$, on obtient
\begin{equation*}
  j = \displaystyle\frac{-(1 + Cr)\ \orr{v}_{ac12} \cdot{} \orr{n}} { 1/m_1 +
    1/m_2 + (\orr{r}_1 \wedge \orr{n})^2 / I_1 + (\orr{r}_2 \wedge \orr{n})^2
    / I_2}
\end{equation*}
et il est possible de déterminer le quadruplet $(\orr{v}'_1, \orr{\omega}'_1,
\orr{v}'_2, \orr{\omega}'_2)$ post-collision. Plus généralement, l'impulsion
résultant d'une collision entre une bactérie $a$ et une bactérie $b$ se notera
$j_{ab}$.

Dans le cas de plusieurs collisions simultanées de la bactérie $1$ avec les
bactéries $(2,\ldots,N)$, dans la mesure ou chaque bactérie est indépendante,
les vitesses linéaires et angulaires post-collision de la bactérie $1$ sont
\begin{equation*}
  \orr{v}'_1 = \orr{v}_{ac1} + \sum^N_{k=2} \frac{j_{1k}\orr{n}}{m_1}
  \qquad
  \orr{\omega}'_1 = \orr{\omega}_{ac1} + \sum^N_{k=2} \frac{\orr{r}_1 \wedge j_{1k}\orr{n}}{I_1}
\end{equation*}
l'impulsion étant le seul terme dépendant de la collision, il est différent pour
chacune des collisions.


% Les places sont limitées
\paragraph{Fonction de transition locale}
La méthode naïve pour détecter les collisions entre un nombre arbitraire de
bactéries est d'effectuer un test de collision pour chaque paire de bactéries,
c'est à dire $\binom{n}{2}$ combinaisons, ce qui n'est pas atteignable pour une
population de l'ordre de \num{E4} individus avec les contraintes que nous nous
sommes fixées. Pour résoudre ce problème, nous utiliserons le fait que, comme
nous effectuons une simulation physique, nous partons du principe que tout ce
qui cause un changement d'état du système a un \emph{effet local} et se
\emph{propage} de proche en proche. C'est bien le chemin que nous avons emprunté
en choisissant les automates cellulaires comme support de la simulation. Nous
devons déterminer comment effectuer notre recherche de collision localement pour
qu'elle soit correcte à l'échelle de la population. Nous allons donc restreindre
la recherche de collision à un petit groupe de bactéries, contenu dans une
cellule de l'automate. Le nombre de bactéries par cellule $\mathit{nbpc}$ dépend
de plusieurs paramètres:
\begin{itemize}
\item l'aire minimum d'une bactérie $S_\mathit{min}$ (dépendant de
  $l_{\mathit{min}}$ et $r_{\mathit{min}}$);
\item l'aire d'une cellule $C_p$ de diffusion du moteur physique.
\end{itemize}
Nous donnons une dimension physique concrète aux cellules car nous allons devoir
superposer les deux automates cellulaires des moteurs physique et chimiques pour
savoir où les bactéries déposent et consomment des morphogènes. De plus, nous
gardons un rapport entier entre les aires des cellules $C_p$ du moteur physique
et $C_c$ du moteur chimique afin que les deux grilles soient toujours alignées.
\[
  \mathit{nbpc} = \left\lfloor
    \displaystyle\frac{C_p}{S_\mathit{min}}
  \right\rfloor
  \qquad \text{et} \qquad
  C_p = n_r \times C_c
\]
Le rapport entre $C_p$ et $C_c$, noté $n_r$, est un entier.

L'espace ainsi segmenté en cellules est une \emph{configuration}. Afin que notre
automate cellulaire soit complet, il nous reste à spécifier la fonction de
transition locale à utiliser pour mettre à jour cette configuration.

\begin{algorithm}
  \SetKwData{Evol}{evol}
  \SetKwData{GCol}{getCollisions}
  \SetKwData{UpBact}{updateBacteria}
  \SetKw{KwFrom}{from}
  \SetKw{And}{and}
  \SetKwData{NBPC}{$\mathit{nbpc}$}
  \SetKwProg{Fn}{Fonction}{ as}

  \Fn{\Evol(c: cell, $n_1$: cell, $\ldots$, $n_8$: cell): cell}{
    \For{$i$ \KwFrom $0$ \KwTo $\NBPC - 1$}{
      \tcp{Calcul de collision}
      $v \gets 0$,
      $\omega \gets 0$\;
      \ForAll{$c_n \in \{c,n_1,\ldots,n_8\}$}{
        \For{$j$ \KwFrom $0$ \KwTo $\NBPC - 1$}{
          \If{$c[i] \ne c_n[j]$ \And $\GCol (c[i],c_n[j],\orr{c},\orr{n},\orr{j})$}{
            $v \gets v + \frac{\orr{j}}{c[i].m}$,
            $\omega \gets \omega + \frac{(\orr{c} - c[i].\orr{p}) \wedge \orr{j}}{c[i].I}$\;
          }
        }
      }

      \tcp{Mise à jour de position}
      $\UpBact (v,\omega)$\;
    }
    \Return{$c'$}\;
  }
  \caption{Fonction de transition locale du moteur physique}
\end{algorithm}

La fonction \ç|getCollisions| effectue les calculs de collision présentés
ci-dessus et retourne \ç|true| si une collision est détectée et \ç|false| sinon.
Les variables \orr{c}, \orr{n} et \orr{j} sont passées par adresse et sont
définies seulement dans le cas d'un collision.

\FloatBarrier
\subsection{Moteur de décision: bactéries et morphogènes réunis}
La coordination entre moteur chimique et moteur physique est le dernier problème
à résoudre. En effet, pour être complet, les deux parties du simulateur doivent
communiquer:
\begin{itemize}
\item les bactéries peuvent \emph{lire} des valeurs des concentrations en
  morphogène \emph{à leur position} et
\item les bactéries peuvent \emph{consommer} ou \emph{déposer} des
  morphogènes à leur position.
\end{itemize}
Lire, consommer et déposer des morphogène implique que la fonction d'évolution
des bactéries doit avoir connaissance de la configuration de l'automate chimique
et doit pouvoir la modifier. Il n'est pas souhaitable que les bactéries écrivent
directement dans les cellules de la configuration de l'automate chimique car ces
valeurs dépendent du temps et sont intégrées numériquement. Une meilleure
solution est de faire déposer par les bactéries une contribution $\epsilon$ qui
sera ajoutée lors de l'intégration numérique pendant la mise à jour de
l'automate chimique. Nous reprenons l'équation \ref{eq:num-int} auquel nous
ajoutons un nouveau terme:
\begin{align*}
  \displaystyle
  \varphi^{t + \Delta t}_{k,i,j} &= \Delta t\,\text{X}_1(\varphi^t_{1,i,j}, \ldots, \varphi^t_{k,i,j})
  + \frac{\text{D}_k\Delta t}{(\Delta x)^2}
  (\varphi^t_{k,i-1,j} + \varphi^t_{k,i+1,j} + \varphi^t_{k,i,j-1} + \varphi^t_{k,i,j+1}
  - 4\varphi^t_{k,i,j}) + \Delta t\,\epsilon_{\varphi_k,i,j}
\end{align*}

Ensuite, les deux configurations des automates physique et chimique doivent se
superposer dans l'espace afin de lier la position d'une bactérie à la
concentration en morphogènes à cette position. À la position du centre d'une
bactérie $(p_x,p_y)$ correspond une cellule de la configuration de l'automate
cellulaire du moteur chimique
\begin{equation*}
  \left(\frac{p_x}{\Delta x},\frac{p_y}{\Delta x}\right)
\end{equation*}
où $\Delta x$ est la taille du côté d'une cellule de cet automate et à condition
qu'à la position $(0,0)$ corresponde la cellule de coordonnées $(0,0)$.

La contrainte suivante concerne la taille des configurations. Nous avons indiqué
dans la section~\ref{sec:boundingbox} que les dimensions des configurations
variaient en fonction de leur contenu. Comme une bactérie doit toujours pouvoir
lire, écrire ou consommer des morphogènes, cela implique que la partie de
l'espace $S_c$ couvert par la configuration du moteur chimique doit inclure la
partie de l'espace $S_p$ couvert par la configuration du moteur physique, soit
$S_p \subset S_c$.


Finalement, lorsque les bactéries ont accès aux concentrations de morphogènes,
elles peuvent s'en servir pour modifier leur comportement ce qui est l'objet du
moteur de décision. Chaque bactérie a un état prédéfini en début de simulation
qui regroupe les différents variables normales et booléennes (dont la valeur par
défaut est donnée entre parenthèses) suivants
\begin{itemize}
\item Division (vrai): si vrai, la bactérie croît et se divise quand sa masse a
  doublé;
\item Décès (faux): si vrai, la bactérie décide d'entrer en apoptose et
  disparaît de la simulation;
\item TumbleRun (vrai): si vrai, la bactérie tourne sur elle-même, sinon la
  bactérie se déplace selon une trajectoire rectiligne uniforme vers l'avant;
\item Production/Consommation du morphogène $\varphi_i$: la bactérie dépose une
  certaine quantité de morphogènes à sa position dans l'environnement. Cette
  quantité est déterminée par le modélisateur en début de simulation.
\end{itemize}
Les transitions entre les états sont également donnés en début de simulation
dans un fichier de description de simulation \sbgp.

Le fonctionnement du moteur de décision est donné par la fonction de transition
décrite dans l'algorithme~\ref{algo:decision}. Cette dernière est appelée sur
chaque cellule de la configuration de l'automate du moteur physique directement,
sans passer par \ppm car il n'y a pas conflit en lecture ou en écriture et
peut être parallélisée naïvement. Cette fonction a de plus à sa disposition
la configuration de l'automate du moteur chimique $\mathit{H}$ ainsi qu'une
configuration $\mathit{E}$ des contributions aux concentrations des morphogènes
de la même taille que $\mathit{H}$. Nous considérerons une simulation où il y
a $N$ morphogènes. Dans l'algorithme~\ref{algo:decision}, toutes les valeurs
précédées de \ç|init| sont issues de l'initialisation de la simulation.


\begin{algorithm}
  \SetKwData{Evol}{evol}
  \SetKw{KwFrom}{from}\SetKw{KwIn}{in}\SetKw{And}{and}\SetKw{KwBreak}{break}
  \SetKwData{NBPC}{$\mathit{nbpc}$}
  \SetKwProg{Fn}{Fonction}{ as}

  \Fn{\Evol(c: cell, H: configuration, E: configuration): cell}{
    \For{$i$ \KwFrom $0$ \KwTo $\NBPC - 1$}{
      \tcp{Indice de la cellule de $H$ correspondant à la bactérie courante}
      $j \gets projection(H,(c[i].px,c[i].py))$\;

      \tcp{Mise à jour de la bactérie pour chaque morphogène}
      \For{$k$ \KwFrom $0$ \KwTo  $N - 1$}{
        E[i][k] <- E[i][k] + \ç|init.behavior.(c[i].type).EmitSignal[k].value|
      }

      \tcp{Transition vers un nouveau type}
      \For{$t$ \KwIn \ç|init.transisition.(c[i].type)|}{
        \If{\ç|init.cond_transition.t|}{
          $c[i].type \gets t$\;
          \KwBreak\;
        }
      }
    }
    \Return{$c'$}\;
  }
  \caption{\label{algo:decision}Fonction de transition locale du moteur de
    décision}
\end{algorithm}


\noindent
En prenant en compte toutes les contraintes précédentes, l'exécution d'un pas de
simulation s'effectue de la façon suivante:
\begin{enumerate}
\item un pas du moteur physique,
\item un pas du moteur chimique,
\item redimensionnement des configurations et
\item un pas du moteur de décision.
\end{enumerate}


\FloatBarrier
\section{Mise en œuvre}\label{sec:mise-en-oeuvre}
Dans cette section nous présentons trois exemples de simulation avec \otb. Pour
chacun de ces exemple nous fournissons le rendu de l'interface graphique et
le fichier SBGP qui l'a engendré. Ces trois exemples mettent en évidence le
fonctionnement de chacun des moteurs (chimique, physique et de décision).

% Chimique BZ
\subsection{Une réaction de Belooussov-Jabotinski}
\label{sec:bz}
\newcommand{\gfp}{\ensuremath{\mathit{bleu}}}%
\newcommand{\rfp}{\ensuremath{\mathit{vert}}}%
Les réactions de Belooussov-Jabotinski (BZ) sont une classe de réactions
chimiques dans laquelle les concentrations de plusieurs réactants oscillent
périodiquement jusqu'à épuisement de l'un d'entre eux. Leur importance réside
dans le fait que l'état d'équilibre n'est pas atteint immédiatement mais après
une série d'oscillations pouvant prendre plusieurs minutes. Quand elle est
effectuée dans une boîte de Petri, il est possible d'observer des vagues de
réactions.

Dans le cas de notre moteur chimique, nous pouvons facilement spécifier les
réactions et les couleurs associés à chacun des morphogènes. Nous construisons
donc cette réaction avec le minimum de réactifs et de réactions nécessaires,
soit deux réactifs \gfp{} et \rfp{} et les trois réactions suivantes:
\begin{align*}
  \gfp + \rfp &\xrightarrow{0.05} 2 \rfp\\
  \gfp        &\xrightarrow{0.05} 2 \gfp\\
  \rfp        &\xrightarrow{0.05} .\\
\end{align*}

\noindent
La configuration de la simulation décrite en \sbgp{} est la suivante
\begin{SBGPcode}
{
  "signals"         : [ ["bleu", 0.0, 1.0], ["vert", 0.0, 1.0] ],
  "reactions"       : [ [["bleu","vert"], ["vert","vert"], 0.05],
                        [["bleu"],        ["bleu","bleu"], 0.05],
                        [["vert"],        [],              0.05] ],
  "type"            : [],
  "behavior"        : {},
  "transition"      : [],
  "cond_transition" : []
}
\end{SBGPcode}
où les quatre derniers champs sont vides car il n'y a pas de bactéries à
simuler. Durant la simulation, il est possible d'observer une alternance entre
les concentrations des deux morphogènes (en vert et en bleu) au cours du temps.
Comme les morphogènes diffusent, ne s'évaporent pas et que l'environnement
est ouvert, les concentrations locales baissent et il devient de plus en
plus difficile de distinguer les deux morphogènes au centre. Par contre, des
vagues de concentrations alternant bleu puis vert se forment à la périphérie
conformément aux réactions BZ effectuées en boîte de Petri.

% Réaction BZ
\begin{figure}[ht]
  \centering
  \includegraphics[width=.245\textwidth]{figures/BZ-reaction1}
  \includegraphics[width=.245\textwidth]{figures/BZ-reaction2}
  \includegraphics[width=.245\textwidth]{figures/BZ-reaction3}
  \includegraphics[width=.245\textwidth]{figures/BZ-reaction4}
  \caption{Quatre états d'une simulation de réaction BZ ordonnés en fonction du
    temps de gauche à droite. Il est possible de noter l'alternance entre le
    vert et le bleu au cours du temps}
  \label{fig:reaction-BZ}
\end{figure}

\subsection{Sectorisation d'une population de \ecoli{}}
Pour la validation du moteur physique, nous avons reproduit une expérience de la
littérature~\cite{hallatschek_genetic_2007}. Dans cet article, les auteurs
étudient la ségrégation génomique à la frontière entre deux populations de
bactéries en croissance.

\paragraph{Protocole expérimental} Les auteurs utilisent deux souches de
bactéries \Ecoli{} de génotypes exactement identique, à l'exception d'une unique
mutation du même gène promoteur induisant
\begin{itemize}
\item la première souche à produire une protéine fluorescente cyan (CFP) et
\item la seconde souche à produire une protéine fluorescente jaune (YFP).
\end{itemize}
On répartit de façon homogène ces deux souches pour obtenir des cultures de
différentes proportions. Une gouttelette de cette culture est posée au centre
d'une gélose pour dénombrement dans un milieu contenant un médium de croissance
enrichi afin d'obtenir une croissance rapide des deux populations. Comme les
bactéries utilisées sont immobiles, le développement de la population s'effectue
aux bords. Dans le cas d'un mélange équilibré contenant 50\% de chaque souche,
il est possible d'observer des secteurs monochromes dans la zone d'expansion de
la population.

D'après les auteurs de~\cite{hallatschek_genetic_2007}, la
\textsc{Fig.}~\ref{fig:sectorisation} à gauche: « […] is a visible manifestation of
random genetic drift acting at the leading edge of a range expansion ». La
dérive génétique (genetic drift) est l'évolution d'une population ou d'une
espèce causée par des phénomènes aléatoires, impossible à prévoir. Du point de
vue génétique, c'est la modification de la fréquence d'un allèle, ou d'un
génotype, au sein d'une population, indépendamment des mutations, de la
sélection naturelle et des migrations.

\paragraph{Simulation avec \otb} Dans notre simulation de cette expérience, nous
avons utilisé deux états différents ayant le même comportement pour les deux
souches de bactéries afin de pouvoir colorer les deux sous populations a
posteriori. Voici le protocole que nous avons suivi:
\begin{enumerate}
\item nous avons initialement laissé croître une population initiale issue d'une
  unique bactérie jusqu'à \num{1054} individus, pour obtenir une goutte de
  l'ordre du dixième de millimètre;
\item nous avons marqué chaque individu par l'état $A$ avec une probabilité de
  50\%, les autres ont été marqué par l'état $B$;
\item Nous avons repris la croissance de cette population jusqu'à atteindre une
  population de \num{271958} individus (voir Fig.~\ref{fig:sectorisation}) après
  vingt minutes de calcul sur un serveur muni d'une carte de calcul
  Tesla\footnote{NVidia Tesla est une carte de calcul dédié munie de plusieurs
    GPU};
\item nous avons fait émettre aux individus de type $A$ un morphogène de couleur
  verte et aux individus type $B$ un morphogène de couleur rouge, afin de
  pouvoir observer qualitativement la répartition des bactéries.
\end{enumerate}

Bien qu'idéalement il nous faudrait une population de l'ordre de \num{E6}
individus pour être dans les mêmes conditions que le protocole expérimental,
nous avons obtenu un résultat qualitativement très proche des photos présentées
dans l'article original. Étonnamment, nos bactéries rudimentaires, par les seules
contraintes mécaniques du moteur physique, se répartissent en secteurs lors de
la croissance de la population.

% Sectorisation
\begin{figure}[ht]
  \centering
  \includegraphics[height=.4\textwidth]{figures/sectors-original}\hspace{1cm}
  \includegraphics[height=.4\textwidth]{figures/sectors}
  \caption[Deux population de bactéries \ecoli{} croissant librement:
    comparaison avec un résultat de la littérature.] {Deux population de
    bactéries \ecoli{} croissant librement dans un milieu riche en nutriment.
    La figure de droite est issue de la fenêtre de rendu du simulateur \otb.
    La figure de gauche est issue de~\cite{hallatschek_genetic_2007} pour
    comparaison}
  \label{fig:sectorisation}
\end{figure}





\FloatBarrier
% Ensemble Cercles concentriques
\subsection{Une population stable ?}
Les bactéries croissent, se divisent et entrent en apoptose. Nous avons tenté
dans cette simulation d'obtenir une population stable ou au moins croissant
linéairement. Intuitivement, l'idée que avons suivie est de faire en sorte que
les bactéries croissent en cercles concentriques autour de la bactérie mère
initiale. Cette dernière dépose un morphogène, appelé «Vert» dans le modèle,
et se divise en deux bactéries: une nouvelle bactérie mère et une bactérie
fille qui change d'état immédiatement. Cette nouvelle bactérie fille dépose
un morphogène «Rouge» dans l'environnement. L'ensemble des bactéries filles
forme alors un premier anneau autour de la bactérie mère. Si ces dernières ne
perçoivent pas un taux de morphogène «Vert» suffisamment grand, signifiant
qu'elles sont trop éloignées, alors elles changent à nouveau d'état.

Voici le comportement que nous avons spécifié en SBGP:
\begin{SBGPcode}
{
  "signals"         : [ ["Vert",  0.01,  2.5],
                        ["Rouge", 0.001, 1.5],
                        ["Bleu",  0.001, 1.5] ],
  "reactions"       : [],
  "type"            : [ "LEADER",
                        "CDAUGHTER",
                        "FDAUGHTER",
                        "DEAD" ],
  "behavior"        : {
    "LEADER" : [ { "Divide" : [] },
                 { "Tumble" : [] },
                 { "EmitSignal" : [ "Vert", 50 ] } ],
    "CDAUGHTER" : [ { "Divide" : 0.005 },
                    { "Tumble"     : [] },
                    { "EmitSignal" : [ "Rouge", 50 ] } ],
    "FDAUGHTER" : [ { "Divide" : 0.005 },
                    { "Tumble"     : [] },
                    { "EmitSignal" : [ "Bleu", 50 ] } ],
    "DEAD" : [ { "Die" : [] } ]
  },
  "transition"      : [
    ["NA","NA","NA","NA"],
    ["CD","NA","NA","NA"],
    ["NA","CF","NA","NA"],
    ["NA","NA","CA","NA"]
  ],
  "cond_transition" : [
    {"CD" : "IsDaughter()"},
    {"CF" : "LessThreshold(Vert,0.1)"},
    {"CA" : "And(LessThreshold(Rouge,0.1),LessThreshold(Bleu,10.0))"}
  ]
}
\end{SBGPcode}

\noindent
Les bactéries peuvent être de quatre types:
\begin{enumerate}
\item \ç|LEADER| est la première bactérie de la simulation. Elle dépose le
  morphogène vert dans l'environnement;
\item \ç|CDAUGHTER| sont les bactéries filles directes de la bactérie de type
  \ç|LEADER|. Elles déposent le morphogène \ç|Rouge| dans l'environnement;
\item \ç|FDAUGHTER| sont les bactéries filles qui sont assez éloignées de la
  bactérie de type \ç|LEADER| d'après la concentration en morphogène \ç|Vert|.
  Elles déposent le morphogène \ç|Bleu| dans l'environnement;
\item \ç|DEAD| sont les bactéries filles qui se sont trop éloignées de la
  bactérie de type \ç|LEADER| et qui meurent quand les concentrations en
  morphogènes \ç|Vert| et \ç|Bleu| sont trop faibles.
\end{enumerate}

\begin{figure}[ht]
  \centering
  \includegraphics{figures/stablePopulation}
  \caption{Graphique de la simulation « Une population stable? » présentant le
    nombre de bactérie en fonction du temps de simulation. Après une phase de
    croissance linéaire, la population se stabilise autour de \num{6000}
    individus}
  \label{graph:stablePopulation}
\end{figure}

\begin{figure}[ht]
  \centering
  \includegraphics[width=.45\textwidth]{figures/ftl01}\hspace{1em}
  \includegraphics[width=.45\textwidth]{figures/ftl02}\vspace{1em}
  \includegraphics[width=.45\textwidth]{figures/ftl03}\hspace{1em}
  \includegraphics[width=.45\textwidth]{figures/ftl05}
  \caption[Différentes étapes de la simulation «Une population
    stable?».]{Différentes étapes de la simulation « Une population stable? ».
    La bactérie mère (en haut à gauche) se divise une première fois et sa fille
    change d'état (en haut à droite). Après un certain nombre de divisions,
    une des bactérie est suffisamment éloignée de la bactérie mère et change
    d'état (en bas à gauche). Plus tard dans la simulation, les bactéries «
    bleues » trop éloignées de la bactérie mère meurent à la périphérie.}
  \label{fig:stablePopulation}
\end{figure}













\FloatBarrier
\section{Conclusion et perspectives}
Cette courte section vise à rappeler les principaux résultats obtenus suite à
nos travaux sur le simulateur de bactéries \otb ainsi que nos futures pistes de
travail sur le sujet.

\subsection{Conclusion}
% Résumé des épisodes précédents
Dans ce chapitre nous présentons \otb, un simulateur d'une population de
bactéries \Ecoli{}. Nous débutons ce chapitre en présentant nos motivations
pour le développement d'\otb ainsi que les choix techniques que nous avons
fait: l'utilisation de \ocaml, de \opencl et \opengl. Nous avons ainsi décidé
de tirer partie du parallélisme présent dans les ordinateurs grand public,
ce qui nous a amené à développer \ppm, un algorithme inspiré des travaux des
chapitres précédents (section~\ref{sec:ppm}). \ppm est conçu de sorte que le
contexte nécessaire à l'exécution de chaque kernel lui soit immédiatement
disponible. Ainsi le travail de chaque kernel est effectué de manière autonome
et en parallèle sans accès concurrent à la mémoire. \ppm est utilisé à la fois
pour l'exécution du moteur chimique et pour l'exécution du moteur physique. Ces
deux moteurs étant indépendants, la conception de \otb repose sur un dernier
moteur de décision permettant aux bactéries et à leur support de communiquer par
le couplage des moteurs chimique et physique. Ces trois moteurs sont présentés
dans la section~\ref{sec:implementation}. Nous clôturons le chapitre par
trois simulations permettant de rendre compte du fonctionnement du simulateur
(section~\ref{sec:mise-en-oeuvre}), montrant le fonctionnement des moteurs
chimique, physique et de décision.

% Résultats importants
\noindent
De nos travaux sur \otb nous tirons deux résultats:
\begin{itemize}
  \item un résultat théorique: l'algorithme \ppm est générique et peut-être
    adapté à d'autres outils travaillant sur les automates cellulaires;
  \item un résultat pratique: le simulateur \otb est un exécutable utilisable
    librement sur différentes plateformes, principalement pour l'étude de la
    morphogénèse dans les populations de \ecoli. Il est modulaire et extensible
    et pourrait donc aisément être adapté à d'autres types de population.
\end{itemize}

\subsection{Perspectives}
Nos travaux futurs s'articulent autour des grands axes suivants, dont
l'objectif principal est la pérennisation d'\otb comme outil de simulation:

\begin{description}

  \item[Calibration]
    Nous projetons de comparer les résultats de simulation avec des résultats
    d'expériences \emph{in vivo} connus afin d'étalonner au mieux les moteurs
    physique et chimique. Si cela s'avère nécessaire, nous changerons
    l'algorithme d'intégration numérique pour une version plus précise, au prix
    de la performance. Nous envisageons d'implémenter bien plus d'exemples de
    modèles qu'actuellement. Nous considérons dans ce but l'établissement d'une
    collaboration avec des scientifiques des sciences du vivant s'intéressant
    à la morphogénèse ou plus largement aux dynamiques d'une population de
    bactéries.

  \item[Performance]
    La performance actuelle du simulateur \otb provient principalement d'une
    planification en amont de sa conception. Il est encore possible d'améliorer
    la vitesse d'exécution du simulateur. Pour cela, il sera nécessaire
    d'effectuer des mesures quantitatives afin d'identifier les parties du
    programme les plus coûteuses en temps de calcul et les optimiser au mieux,
    par rapport à l'environnement dans lequel le simulateur est utilisé et
    en fonction des ressources en espace disponible sur le périphérique de
    calcul. Il pourra ensuite être envisageable de comparer les performances de
    \otb par rapport à d'autres simulateurs de population de bactéries comme
    \gro~\cite{jang_specification_2012}.

  \item[Interface de sortie]
    L'interface utilisateur de \otb est à ce jour rudimentaire. Pour qu'un plus
    grand nombre d'utilisateurs puissent utiliser le simulateur au quotidien
    il est impératif de lui ajouter quelques fonctionnalités: un meilleur
    contrôle à la souris, l'ajout d'obstacles et de flux dans l'environnement
    afin de tester l'influence de contraintes mécaniques sur une population
    en croissance ou encore présenter plusieurs vues en simultané de la même
    population. Il est envisageable d'ajouter d'autres modes de sortie, comme
    des mesures quantitatives agrégées sous forme de graphes en fonction du
    temps de la simulation comme, par exemple, l'orientation moyenne des
    bactéries, leur nombre, leur taille moyenne, leur volume moyen, etc.

  \item[Interface d'entrée]
    La description d'une simulation par un fichier \sbgp{} peut encore être
    améliorée. Il est d'abord nécessaire d'optimiser le contenu de ce fichier
    en fonction des informations réellement nécessaire à une simulation \otb
    ce qui sera possible une fois le point «Calibration» effectué et après
    discussion avec les futurs utilisateurs d'\otb. Il est aussi envisageable
    de contrôler \otb à partir d'un programme externe, par exemple une
    interface graphique regroupant les interfaces d'entrée et de sortie pour un
    usage plus intuitif.

  \item[Formes de bactérie]
    De part sa conception, \otb est modulaire et extensible. Une piste
    de travail futur est l'exploration de la simulation d'autre formes de
    bactéries du genre bacillus (comme \ecoli), diplobacillus (des couples de
    bâtonnets), streptobacillus (des chaînes de bacillus), coccus (des sphères),
    coccobacillus (des ovales), etc. Changer de forme de bactérie implique de
    mettre à jour l'algorithme de collision qui a été développé spécifiquement
    pour les bactéries bacillaires. Ces développements se feront également en
    lien avec le point «Calibration» en fonction des besoins.

  \item[Grille hexagonale]
    Une propriété des grilles hexagonales est que toutes les cellules ont un
    bord en commun et leur centres sont équidistants (ce qui n'est pas le cas
    sur une grille carrée, les centres des cellules en diagonale sont plus
    éloignés). Nous souhaitons explorer le fonctionnement de \ppm sur ce support
    afin de déterminer s'il est possible d'obtenir une diminution des erreurs
    de la méthode d'intégration actuelle du moteur chimique et du calcul de
    collision du moteur physique.

  \item[\mgs]
    Dans le point «Interface d'entrée» nous avons indiqué qu'il était
    envisageable qu'un programme externe serve d'interface de contrôle pour
    \otb. \mgs est un bon candidat pour cette tâche. En effet, nous pouvons
    considérer que \otb définit une collection topologique dont les éléments
    sont les bactéries. Cette collection est paramétrée par les réactions
    entre les morphogènes. Chaque élément de cette collection a pour valeur
    la spécification d'un réseau de régulation génétique sous la forme d'un
    automate avec nombre fini d'états.

\end{description}