\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}